Generating final model predictions after relativizing

In this document, we integrate our forest classification model as well as models for each different functional group and ecoregion we considered into CONUS-wide predictions of both absolute and relative cover. We show maps of these model predictions that are generated using contemporary climate, weather, and soils data, and assess how these predictions compare to observed data using maps of residuals, quantile plots, bias, and RMSE. We also show model predictions of both absolute and relative cover made using modeled climate and weather data from the end of the 20th century under RCP 8.5. We used data from two different generalized circulation models (GCMs) from CMIP5: “BNU-ESM”, a cooler and wetter model, and “IPSL-CM5A-MR (France)”, a hotter and drier model. For these future predictions of cover, we display maps that show how these future model predictions compare to modeled contemporary cover. We present all of this information (maps of contemporary cover predictions, residuals, quantile plots, RMSE, bias, and maps of future cover based on GCMs) for both absolute and relative cover.

Setup the environment:

library(tidyverse)
library(sf)
library(terra)
library(kableExtra)
library(knitr)
library(USA.state.boundaries)
library(tidyterra)
library(ggpubr)
library(USA.state.boundaries)
library(ggrepel)
# set ggplot theme default 
theme_set(theme_classic())

# Load User defined parameters and functions
print(params)

## $readParams
## [1] TRUE

# set to true if want to run for a limited number of rows (i.e. for code testing)
readParams <- params$readParams
## function to maek predictionss f
makePredictions <- function(predictionDF, modelObject, scale  = TRUE) {
  ##
  # predictionDF <- climDatPred
  # modelObject <- bestLambdaMod_grassShrub_totalHerb
  # # ##
  # predictionDF = modDat_quantile,
  # modelObject = forb_CONUS, scale = FALSE
  # 
  # pretend to scale the input variables so the model object can predict accurately
  if(scale == TRUE) {
  predictionDF <- predictionDF %>% 
  mutate(across(all_of(prednames), base::scale,scale = FALSE, center = FALSE)) 
  }
  # modelPredictions
  modelPreds <- predict(object = modelObject, newdata = predictionDF, type = "response")
  # add predictions back into the input data.frame
  predictionDF <- predictionDF %>% 
    cbind(modelPreds)
  
  # truncate all predictions to max out at 100 
  #predictionDF[predictionDF$modelPreds>100 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 100
predictionDF[predictionDF$modelPreds < 0 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 0

  # print predicted data
 return(predictionDF)
}

## source functions for quantile plots
source("../../../Functions/glmTransformsIterates.R")
source("../../../Functions/transformPreds.R")

# Load Data ---------------------------------------------------------------
# data ready for cover modeling (that has been scaled)
modDat_1_s <- readRDS("./models/scaledModelInputData.rds")
# data ready for yes/no forest modeling (that has been scaled)
modDat_tree_s <- readRDS("./models/scaledModelInputData_yesNoTreeMod.rds")
  
# get the soil raster, which we'll use for rasterizing the imagery
soilRastTemp <- readRDS("../../../Data_processed/SoilsRaster.rds") %>% 
terra::unwrap()

# make a map of the predictions
test_rast <-  rast("../../../Data_raw/dayMet/rawMonthlyData/orders/70e0da02b9d2d6e8faa8c97d211f3546/Daymet_Monthly_V4R1/data/daymet_v4_prcp_monttl_na_1980.tif") %>% 
  terra::aggregate(fact = 8, fun = "mean", na.rm = TRUE)

## |---------|---------|---------|---------|=========================================                                          

# download map info for visualization
data(state_boundaries_wgs84) 

cropped_states <- suppressMessages(state_boundaries_wgs84 %>%
  dplyr::filter(NAME!="Hawaii") %>%
  dplyr::filter(NAME!="Alaska") %>%
  dplyr::filter(NAME!="Puerto Rico") %>%
  dplyr::filter(NAME!="American Samoa") %>%
  dplyr::filter(NAME!="Guam") %>%
  dplyr::filter(NAME!="Commonwealth of the Northern Mariana Islands") %>%
  dplyr::filter(NAME!="United States Virgin Islands") %>%
  dplyr::filter(NAME!= "U.S. Virgin Islands") %>% 
  sf::st_sf() %>%
  sf::st_transform(sf::st_crs(test_rast))) 
  #sf::st_crop(sf::st_bbox(modDat_1_sf)+c(-1,-1,1,1))

## add ecoregion boundaries (for our ecoregion level model)
regions <- sf::st_read(dsn = "../../../Data_raw/Level2Ecoregions/", layer = "NA_CEC_Eco_Level2") 

## Reading layer `NA_CEC_Eco_Level2' from data source 
##   `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/cleanPED/PED_vegClimModels/Data_raw/Level2Ecoregions' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area

regions <- regions %>% 
  st_transform(crs = st_crs(test_rast)) %>% 
  st_make_valid() 
ecoregionLU <- data.frame("NA_L1NAME" = sort(unique(regions$NA_L1NAME)), 
                        "newRegion" = c(NA, "Forest", "dryShrubGrass", 
                                        "dryShrubGrass", "Forest", "dryShrubGrass",
                                       "dryShrubGrass", "Forest", "Forest", 
                                       "dryShrubGrass", "Forest", "Forest", 
                                       "Forest", "Forest", "dryShrubGrass", 
                                       NA
                                        ))
goodRegions <- regions %>% 
  left_join(ecoregionLU)
mapRegions <- goodRegions %>% 
  filter(!is.na(newRegion)) %>% 
  group_by(newRegion) %>% 
  summarise(geometry = sf::st_union(geometry)) %>% 
  ungroup() %>% 
  st_simplify(dTolerance = 1000)

## contemporary climate data
climDat_temp <- readRDS("../../../Data_processed/EcoRegion_climSoilData.rds")

climDat_timeTemp <- climDat_temp %>% 
  mutate(spatialID = paste0(round(x,4), "_", round(y,4))) %>% 
  drop_na()

spatID_lu <- data.frame("spatialID" = unique(climDat_timeTemp$spatialID), 
                        "uniqueID" = 1:length(unique(climDat_timeTemp$spatialID)))

climDat_timeTemp <- climDat_timeTemp %>% 
  left_join(spatID_lu)

## because of the spatial joining we did, there are some "locations" that have data for less than the true time span we want (2010-2020), or multiple observations per year... remove those so we have the even representation over the time span we're looking for 
set.seed("12011993")
goodTable <- data.frame(table(climDat_timeTemp$uniqueID))
goodIDs_temp <- as.numeric(goodTable[goodTable$Freq == 11,"Var1"])

#goodIDs <- unique(names(table(climDat_timeTemp$uniqueID))[(table(climDat_timeTemp$uniqueID)==11)])
# select 1000 random IDs
goodIDs <- goodIDs_temp[sample(x = 1:length(goodIDs_temp), size = 1000, replace = FALSE)]

## get climate data for 500 different points randomly chosen across CONUS for all years in the contemporary climate/weather/soils dataset 
climDat_time <- climDat_timeTemp %>% 
  filter(uniqueID %in% goodIDs) %>% 
  dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE, 
         NA_L1NAME, NA_L1KEY, year, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity, uniqueID, spatialID) %>% 
rename("tmin" = tmin_meanAnnAvg_CLIM, 
     "tmax" = tmax_meanAnnAvg_CLIM, #1 
     "tmean" = tmean_meanAnnAvg_CLIM, 
     "prcp" = prcp_meanAnnTotal_CLIM, 
     "t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
     "t_cold" = T_coldestMonth_meanAnnAvg_CLIM, 
     "prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
     "prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM, 
     "prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
     "prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM,  #3
     "abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM, 
     "isothermality" = isothermality_meanAnnAvg_CLIM, #4
     "annWatDef" = annWaterDeficit_meanAnnAvg_CLIM, 
     "annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
     "VPD_mean" = annVPD_mean_meanAnnAvg_CLIM, 
     "VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
     "VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
     "VPD_max_95" = annVPD_max_95percentile_CLIM, 
     "annWatDef_95" = annWaterDeficit_95percentile_CLIM, 
     "annWetDegDays_5" = annWetDegDays_5percentile_CLIM, 
     "frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM, 
     "frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM, 
     "soilDepth" = soilDepth, #7
     "clay" = surfaceClay_perc, 
     "sand" = avgSandPerc_acrossDepth, #8
     "coarse" = avgCoarsePerc_acrossDepth, #9
     "carbon" = avgOrganicCarbonPerc_0_3cm, #10
     "AWHC" = totalAvailableWaterHoldingCapacity,
     ## anomaly variables
     tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
     tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
     tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
    prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
      t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom,  #19
     t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
      prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
      precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom,  #22
    prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23 
     prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
      aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25  
    isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
       annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
     annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom,  #28
      VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
      VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom,  #30
      VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom,  #31
     VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
      annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33 
      annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom ,  #34
    frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35 
      frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
  ) %>% 
  dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr)) 


# rename
climDat <- climDat_temp %>% 
  filter(year == 2016) %>% 
  dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE, 
         NA_L1NAME, NA_L1KEY, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity) %>% 
rename("tmin" = tmin_meanAnnAvg_CLIM, 
     "tmax" = tmax_meanAnnAvg_CLIM, #1 
     "tmean" = tmean_meanAnnAvg_CLIM, 
     "prcp" = prcp_meanAnnTotal_CLIM, 
     "t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
     "t_cold" = T_coldestMonth_meanAnnAvg_CLIM, 
     "prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
     "prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM, 
     "prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
     "prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM,  #3
     "abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM, 
     "isothermality" = isothermality_meanAnnAvg_CLIM, #4
     "annWatDef" = annWaterDeficit_meanAnnAvg_CLIM, 
     "annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
     "VPD_mean" = annVPD_mean_meanAnnAvg_CLIM, 
     "VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
     "VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
     "VPD_max_95" = annVPD_max_95percentile_CLIM, 
     "annWatDef_95" = annWaterDeficit_95percentile_CLIM, 
     "annWetDegDays_5" = annWetDegDays_5percentile_CLIM, 
     "frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM, 
     "frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM, 
     "soilDepth" = soilDepth, #7
     "clay" = surfaceClay_perc, 
     "sand" = avgSandPerc_acrossDepth, #8
     "coarse" = avgCoarsePerc_acrossDepth, #9
     "carbon" = avgOrganicCarbonPerc_0_3cm, #10
     "AWHC" = totalAvailableWaterHoldingCapacity,
     ## anomaly variables
     tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
     tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
     tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
    prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
      t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom,  #19
     t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
      prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
      precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom,  #22
    prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23 
     prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
      aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25  
    isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
       annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
     annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom,  #28
      VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
      VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom,  #30
      VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom,  #31
     VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
      annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33 
      annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom ,  #34
    frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35 
      frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
  ) %>% 
  dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr))

rm(climDat_temp) 
gc()

##              used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
## Ncells    3109203  166.1    5111205   273.0         NA    3965807   211.8
## Vcells 1074653381 8199.0 3081476488 23509.8      65536 2842501221 21686.6

## now, scale the contemporary climate data for use in cover models 
# get the scaling factors 
scaleParams <- modDat_1_s %>% 
  #filter(Year == 2016) %>% 
  dplyr::select(tmin_s:AWHC_s) %>% 
  reframe(across(all_of(names(.)), attributes)) 

# apply the scaling factors to the contemporary climate data 
namesToScale <- climDat %>% 
  dplyr::select(tmin:frostFreeDays, tmean_anom:frostFreeDays_anom, soilDepth:AWHC) %>% 
  names()

climDat_scaled <- map(namesToScale, .f = function(x) {
  x_new <- (climDat[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")

climDatPred <- climDat %>% 
  dplyr::select(NA_L1CODE:y) %>% 
  cbind(climDat_scaled)
names(climDatPred)[7:56] <- str_remove(names(climDatPred)[7:56], pattern = "_s$")

# apply the scaling factors to the contemporary climate data used for predictions over time 
climDat_scaled <- map(namesToScale, .f = function(x) {
  x_new <- (climDat_time[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")

climDat_time_Pred <- climDat_time %>% 
  dplyr::select(NA_L1CODE:y,spatialID) %>% 
  cbind(climDat_scaled)
names(climDat_time_Pred)[9:58] <- str_remove(names(climDat_time_Pred)[9:58], pattern = "_s$")

rm(climDat_scaled) 
gc()

##              used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
## Ncells    3116902  166.5    5111205   273.0         NA    5111205   273.0
## Vcells 1100219206 8394.1 3081476488 23509.8      65536 2842501221 21686.6

prednames_s <-  modDat_1_s %>%
  dplyr::select(tmin_s:AWHC_s) %>%
  names()
prednames <- str_replace(prednames_s, pattern = "_s$", replacement = "")

## read in scaled data for future climate 
 forecastClimSoilsDatPred_1 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5.rds")
 forecastClimSoilsDatPred_2 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5.rds")
  # read in unscaled data for future climate
 forecastClimSoilsDat_1 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5_UNSCALED.rds")
forecastClimSoilsDat_2 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5_UNSCALED.rds")

# now, use the scaling factors we used for the forest model to scale the contemporary and forecasted climate data to predict with that model
## now, scale the contemporary climate data for use in cover models 
# get the scaling factors 
scaleParams_tree <- modDat_tree_s %>% 
  #filter(Year == 2016) %>% 
  dplyr::select(tmin_s:AWHC_s) %>% 
  reframe(across(all_of(names(.)), attributes)) 

# apply the scaling factors to the contemporary climate data 
namesToScale_tree <- climDat %>% 
  dplyr::select(tmin:frostFreeDays, tmean_anom:frostFreeDays_anom, soilDepth:AWHC) %>% 
  names()

climDat_scaled_tree <- map(namesToScale_tree, .f = function(x) {
  x_new <- (climDat[,x] - scaleParams_tree[,paste0(x, "_s")]$`scaled:center`)/scaleParams_tree[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(climDat_scaled_tree) <- paste0(namesToScale_tree, "_s")

climDatPred_tree <- climDat %>% 
  dplyr::select(NA_L1CODE:y) %>% 
  cbind(climDat_scaled_tree)
names(climDatPred_tree)[7:56] <- str_remove(names(climDatPred_tree)[7:56], pattern = "_s$")

# apply the scaling factors to the contemporary climate data used for predictions over time 
climDat_time_scaled_tree <- map(namesToScale_tree, .f = function(x) {
  x_new <- (climDat_time[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(climDat_time_scaled_tree) <- paste0(namesToScale_tree, "_s")

climDat_time_Pred_tree <- climDat_time %>% 
  dplyr::select(NA_L1CODE:y,spatialID) %>% 
  cbind(climDat_time_scaled_tree)
names(climDat_time_Pred_tree)[9:58] <- str_remove(names(climDat_time_Pred_tree)[9:58], pattern = "_s$")

rm(climDat_scaled_tree, climDat_time_scaled_tree) 
gc()

##              used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
## Ncells    3117358  166.5    5111205   273.0         NA    5111205   273.0
## Vcells 1136571051 8671.4 3081476488 23509.8      65536 2842501221 21686.6

# now scale the future climate data
# forecastClimSoilsDat_1
forecastClimSoilsDat_scaledTree_1 <- map(namesToScale_tree, .f = function(x) {
  x_new <- (forecastClimSoilsDat_1[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(forecastClimSoilsDat_scaledTree_1) <- paste0(namesToScale_tree, "_s")

forecastClimSoilsDat_scaledTree_1 <- forecastClimSoilsDat_1 %>% 
  dplyr::select(x:Year,newRegion) %>% 
  cbind(forecastClimSoilsDat_scaledTree_1)
names(forecastClimSoilsDat_scaledTree_1)[5:54] <- str_remove(names(forecastClimSoilsDat_scaledTree_1)[5:54], pattern = "_s$")

# forecastClimSoilsDat_2
forecastClimSoilsDat_scaledTree_2 <- map(namesToScale_tree, .f = function(x) {
  x_new <- (forecastClimSoilsDat_2[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(forecastClimSoilsDat_scaledTree_2) <- paste0(namesToScale_tree, "_s")

forecastClimSoilsDat_scaledTree_2 <- forecastClimSoilsDat_2 %>% 
  dplyr::select(x:Year,newRegion) %>% 
  cbind(forecastClimSoilsDat_scaledTree_2)
names(forecastClimSoilsDat_scaledTree_2)[5:54] <- str_remove(names(forecastClimSoilsDat_scaledTree_2)[5:54], pattern = "_s$")

First, we show predictions from the model of forest classification

These maps show predictions of forest classification using contemporary and modeled future climate data, where the predicted value is the forest classification. We arrive at these classification values by translating the continuous probability of forest that was predicted by a binomial generalized linear model to binary values using a cutoff of 0.382 (values greater than 0.382 are forest, while values less than or equal to 0.382 are not forest). This threshold is derived from the average prevalence of forest in the data used to fit the binomial model. Note that while the climate data used in model predictions is either contemporary or from GCMs, the soils data used is static across time.

forestMod <- fit_glm_bestLambda <- readRDS("./models/yesOrNoTrees_bestLambdaGLM_10.rds")

## print out model statement
coefficientTable <- data.frame(coefficients(forestMod))
responseVar <- "P(forest)"
coefficientTable$coefficientName <- rownames(coefficientTable)
coefficientTable$coefficients.forestMod. <- round(coefficientTable$coefficients.forestMod., 4)
  
  # print out coefficients w/ coefficient names
tempNames <- paste0(
  apply(coefficientTable, MARGIN = 1, FUN = function(x) {
    if (x["coefficientName"] == "(Intercept)") {
      paste0(x["coefficients.forestMod."])
    } else {  
      paste0(x["coefficients.forestMod."], "*", x["coefficientName"])
    }
    }
  ),
  collapse = " + "
)

# print the unscaled model statement
  unscaledModelName <- paste0(responseVar, "~ 1/ (1 + exp(-(", tempNames,")))")
 #print(unscaledModelName)

# Here are maps of the ecoregion model predictions with both contemporary and modeled future climate data 

## get model data used for ecoregion fitting
# modDat_ecoregionFit <- modDat_tree_s
# 
# modDat_ecoregionFit <- modDat_ecoregionFit %>% 
#   rename("Long" = x, "Lat" = y) %>% 
#   dplyr::select(c(newRegion, #swe_meanAnnAvg_CLIM:
#                   tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_CLIM,
#                                          soilDepth                  , surfaceClay_perc          ,                
#                                          avgSandPerc_acrossDepth    ,  avgCoarsePerc_acrossDepth,                 
#                                          avgOrganicCarbonPerc_0_3cm ,  totalAvailableWaterHoldingCapacity,
#                                          Long, Lat)) %>% 
#   mutate(newRegion = as.factor(newRegion)) %>% 
#   drop_na()

# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDatPred_unscaled <- climDat

#rm(climDat) 
#gc()
preds_byHand <- climDatPred_unscaled %>% 
  mutate(pred =  as.numeric(1/ (1 + exp(-(-0.67478165 +  
                                            1.71345065*((prcp - 483.34337211) / 321.52170566) + 
                                            -0.01633932*((prcp_seasonality - 0.97338094) / 0.23365057) + 
                                            -0.91285243*((prcpTempCorr - 0.01681097) / 0.39882063) +  
                                            0.53377439*((annWetDegDays - 1808.61043539) / 1083.25298131) + 
                                            -0.60523936*((annWatDef_95 - 164.68611646) / 104.14892564) +  
                                            0.20578879*((coarse - 9.66120017) / 10.16883531) + 
                                            -0.61568228*((AWHC - 14.61190409) / 5.65765542) + 
                                            -0.31818631*I(((tmean - 11.0291275) / 4.89048541)^2) + 
                                            0.20104126*I(((prcp_seasonality - 0.97338094) / 0.23365057)^2) + 
                                            -0.89055646*I(((prcpTempCorr - 0.01681097) / 0.39882063)^2) +  
                                            0.14226917*I(((isothermality - 37.76731082) / 5.2256935)^2) + 
                                            -0.11810211*I(((clay - 19.59171509) / 8.51944876)^2) +  
                                            0.11278028*I(((AWHC - 14.61190409) / 5.65765542)^2) +  
                                            0.91656456*((annWetDegDays - 1808.61043539) / 1083.25298131)*((annWatDef_95 - 164.68611646) / 104.14892564) + 
                                            -0.29650844*((prcp_dry - 3.18622575) / 5.36400552)*((isothermality - 37.76731082) / 5.2256935) + 
                                            -0.02061048*((prcpTempCorr - 0.01681097) / 0.39882063)*((isothermality - 37.76731082) / 5.2256935) +
                                            -0.28504255*((isothermality - 37.76731082) / 5.2256935)*((tmean - 11.0291275) / 4.89048541) + 
                                            0.50549093*((prcpTempCorr - 0.01681097) / 0.39882063)*((prcp_dry - 3.18622575) / 5.36400552) +
                                            -0.25285633*((prcpTempCorr - 0.01681097) / 0.39882063)*((tmean - 11.0291275) / 4.89048541) + 
                                            -0.16290284*((AWHC - 14.61190409) / 5.65765542)*((clay - 19.59171509) / 8.51944876) +  
                                            0.09163249*((coarse - 9.66120017) / 10.16883531)*((clay - 19.59171509) / 8.51944876))))))
# convert the continuous probability predictions to binary using previously-identified threshold 
preds_byHand <- preds_byHand %>% 
  mutate(pred_binary = pred,
    pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
         pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))
# ggplot(preds_byHand) +
#   geom_point(aes(x = x, y = y, col = pred_binary))
#predictionsModel <- predict(object = forestMod, type = "response", newdata = climDatPred_unscaled)
# 
# 
# ## prediction w/ data used to fit
# preds_auto_rawDat <- modDat_tree_s %>% 
#   select(x, y, Year, prcp_s, prcp_seasonality_s, prcpTempCorr_s, annWetDegDays_s, annWatDef_95_s, coarse_s, AWHC_s, tmean_s, isothermality_s, prcp_dry_s, clay_s) %>% 
#   rename(prcp = prcp_s, 
#          prcp_seasonality = prcp_seasonality_s, 
#         prcpTempCorr = prcpTempCorr_s, 
#         annWetDegDays = annWetDegDays_s, 
#         annWatDef_95 = annWatDef_95_s, 
#         coarse = coarse_s, 
#         AWHC = AWHC_s,
#         tmean = tmean_s, 
#         isothermality = isothermality_s, 
#         prcp_dry = prcp_dry_s, 
#         clay = clay_s
#         ) 
#  preds_auto_rawDat$pred <- predict(object = forestMod, newdata = preds_auto_rawDat, type = "response")
#  preds_auto_rawDat <- preds_auto_rawDat %>% 
#   mutate(pred_binary = pred,
#     pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
#          pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))
# 
# # plot(density(-0.25285633*((modDat_tree_s$prcpTempCorr - 0.01681097) / 0.39882063)*((modDat_tree_s$tmean - 11.0291275) / 4.89048541) ))
# # lines(density(-0.25285633*(((modDat_tree_s$prcpTempCorr - 0.01681097) / 0.39882063)*((modDat_tree_s$tmean - 11.0291275) / 4.89048541)) ), col = "red")
# # 
# # lines(density((modDat_tree_s[!is.na(modDat_tree_s$clay),]$clay - 19.59171509) / 8.51944876), col = "red")
# 
# preds_byHand_rawDat <- modDat_tree_s %>% 
#   drop_na(prcp, prcp_seasonality, prcpTempCorr, annWetDegDays, annWatDef_95, coarse, AWHC, tmean, isothermality, prcp_dry, clay) %>% 
#   mutate(pred =  as.numeric(1/ (1 + exp(-(-0.67478165 +  
#                                             1.71345065*((prcp - 483.34337211) / 321.52170566) + 
#                                             -0.01633932*((prcp_seasonality - 0.97338094) / 0.23365057) + 
#                                             -0.91285243*((prcpTempCorr - 0.01681097) / 0.39882063) +  
#                                             0.53377439*((annWetDegDays - 1808.61043539) / 1083.25298131) + 
#                                             -0.60523936*((annWatDef_95 - 164.68611646) / 104.14892564) +  
#                                             0.20578879*((coarse - 9.66120017) / 10.16883531) + 
#                                             -0.61568228*((AWHC - 14.61190409) / 5.65765542) + 
#                                             -0.31818631*I(((tmean - 11.0291275) / 4.89048541)^2) + 
#                                             0.20104126*I(((prcp_seasonality - 0.97338094) / 0.23365057)^2) + 
#                                             -0.89055646*I(((prcpTempCorr - 0.01681097) / 0.39882063)^2) +  
#                                             0.14226917*I(((isothermality - 37.76731082) / 5.2256935)^2) + 
#                                             -0.11810211*I(((clay - 19.59171509) / 8.51944876)^2) +  
#                                             0.11278028*I(((AWHC - 14.61190409) / 5.65765542)^2) +  
#                                             0.91656456*((annWetDegDays - 1808.61043539) / 1083.25298131)*((annWatDef_95 - 164.68611646) / 104.14892564) + 
#                                             -0.29650844*((prcp_dry - 3.18622575) / 5.36400552)*((isothermality - 37.76731082) / 5.2256935) + 
#                                             -0.02061048*((prcpTempCorr - 0.01681097) / 0.39882063)*((isothermality - 37.76731082) / 5.2256935) +
#                                             -0.28504255*((isothermality - 37.76731082) / 5.2256935)*((tmean - 11.0291275) / 4.89048541) + 
#                                             0.50549093*((prcpTempCorr - 0.01681097) / 0.39882063)*((prcp_dry - 3.18622575) / 5.36400552) +
#                                             -0.25285633*((prcpTempCorr - 0.01681097) / 0.39882063)*((tmean - 11.0291275) / 4.89048541) + 
#                                             -0.16290284*((AWHC - 14.61190409) / 5.65765542)*((clay - 19.59171509) / 8.51944876) +  
#                                             0.09163249*((coarse - 9.66120017) / 10.16883531)*((clay - 19.59171509) / 8.51944876)
#                                           )
#                                         )
#                                 )))
# # convert the continuous probability predictions to binary using previously-identified threshold 
# preds_byHand_rawDat <- preds_byHand_rawDat %>% 
#   mutate(pred_binary = pred,
#     pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
#          pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))
# 

# predict with modeled future climate data v1
preds_future1_byHand <- forecastClimSoilsDat_1 %>% 
  mutate(pred =  as.numeric(1/ (1 + exp(-(-0.67478165 +  
                                            1.71345065*((prcp - 483.34337211) / 321.52170566) + 
                                            -0.01633932*((prcp_seasonality - 0.97338094) / 0.23365057) + 
                                            -0.91285243*((prcpTempCorr - 0.01681097) / 0.39882063) +  
                                            0.53377439*((annWetDegDays - 1808.61043539) / 1083.25298131) + 
                                            -0.60523936*((annWatDef_95 - 164.68611646) / 104.14892564) +  
                                            0.20578879*((coarse - 9.66120017) / 10.16883531) + 
                                            -0.61568228*((AWHC - 14.61190409) / 5.65765542) + 
                                            -0.31818631*I(((tmean - 11.0291275) / 4.89048541)^2) + 
                                            0.20104126*I(((prcp_seasonality - 0.97338094) / 0.23365057)^2) + 
                                            -0.89055646*I(((prcpTempCorr - 0.01681097) / 0.39882063)^2) +  
                                            0.14226917*I(((isothermality - 37.76731082) / 5.2256935)^2) + 
                                            -0.11810211*I(((clay - 19.59171509) / 8.51944876)^2) +  
                                            0.11278028*I(((AWHC - 14.61190409) / 5.65765542)^2) +  
                                            0.91656456*((annWetDegDays - 1808.61043539) / 1083.25298131)*((annWatDef_95 - 164.68611646) / 104.14892564) + 
                                            -0.29650844*((prcp_dry - 3.18622575) / 5.36400552)*((isothermality - 37.76731082) / 5.2256935) + 
                                            -0.02061048*((prcpTempCorr - 0.01681097) / 0.39882063)*((isothermality - 37.76731082) / 5.2256935) +
                                            -0.28504255*((isothermality - 37.76731082) / 5.2256935)*((tmean - 11.0291275) / 4.89048541) + 
                                            0.50549093*((prcpTempCorr - 0.01681097) / 0.39882063)*((prcp_dry - 3.18622575) / 5.36400552) +
                                            -0.25285633*((prcpTempCorr - 0.01681097) / 0.39882063)*((tmean - 11.0291275) / 4.89048541) + 
                                            -0.16290284*((AWHC - 14.61190409) / 5.65765542)*((clay - 19.59171509) / 8.51944876) +  
                                            0.09163249*((coarse - 9.66120017) / 10.16883531)*((clay - 19.59171509) / 8.51944876))))))
preds_future1_byHand <- preds_future1_byHand %>% 
  mutate(pred_binary = pred,
    pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
         pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))
# predict with modeled future climate data v2
preds_future2_byHand <- forecastClimSoilsDat_2 %>% 
  mutate(pred =  as.numeric(1/ (1 + exp(-(-0.67478165 +  
                                            1.71345065*((prcp - 483.34337211) / 321.52170566) + 
                                            -0.01633932*((prcp_seasonality - 0.97338094) / 0.23365057) + 
                                            -0.91285243*((prcpTempCorr - 0.01681097) / 0.39882063) +  
                                            0.53377439*((annWetDegDays - 1808.61043539) / 1083.25298131) + 
                                            -0.60523936*((annWatDef_95 - 164.68611646) / 104.14892564) +  
                                            0.20578879*((coarse - 9.66120017) / 10.16883531) + 
                                            -0.61568228*((AWHC - 14.61190409) / 5.65765542) + 
                                            -0.31818631*I(((tmean - 11.0291275) / 4.89048541)^2) + 
                                            0.20104126*I(((prcp_seasonality - 0.97338094) / 0.23365057)^2) + 
                                            -0.89055646*I(((prcpTempCorr - 0.01681097) / 0.39882063)^2) +  
                                            0.14226917*I(((isothermality - 37.76731082) / 5.2256935)^2) + 
                                            -0.11810211*I(((clay - 19.59171509) / 8.51944876)^2) +  
                                            0.11278028*I(((AWHC - 14.61190409) / 5.65765542)^2) +  
                                            0.91656456*((annWetDegDays - 1808.61043539) / 1083.25298131)*((annWatDef_95 - 164.68611646) / 104.14892564) + 
                                            -0.29650844*((prcp_dry - 3.18622575) / 5.36400552)*((isothermality - 37.76731082) / 5.2256935) + 
                                            -0.02061048*((prcpTempCorr - 0.01681097) / 0.39882063)*((isothermality - 37.76731082) / 5.2256935) +
                                            -0.28504255*((isothermality - 37.76731082) / 5.2256935)*((tmean - 11.0291275) / 4.89048541) + 
                                            0.50549093*((prcpTempCorr - 0.01681097) / 0.39882063)*((prcp_dry - 3.18622575) / 5.36400552) +
                                            -0.25285633*((prcpTempCorr - 0.01681097) / 0.39882063)*((tmean - 11.0291275) / 4.89048541) + 
                                            -0.16290284*((AWHC - 14.61190409) / 5.65765542)*((clay - 19.59171509) / 8.51944876) +  
                                            0.09163249*((coarse - 9.66120017) / 10.16883531)*((clay - 19.59171509) / 8.51944876)
                                          )))))
preds_future2_byHand <- preds_future2_byHand %>% 
  mutate(pred_binary = pred,
    pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
         pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))

## Plot predictions with contemporary data
# make into a raster
plotObs <- preds_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "pred_binary", 
                   fun = mean, na.rm = TRUE) 

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)

plotObs_2 <- plotObs %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_current <- ggplot() +
  geom_spatraster(data = plotObs_2, maxcell = 50000000) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = c("Forest classification model predictions using 
                    contemporary dayMet data (from 2016)"),
     subtitle = c("Values are the average binary predictions within each 8km x 8 km grid cell")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "", "", "", "forest")) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

# # for predictions w/ raw data
# plotObs_raw <- preds_byHand_rawDat %>%
#          #drop_na(paste(response)) %>%
#   st_as_sf(coords = c("x", "y"), crs = crs("EPSG:4326")) %>%
#   st_transform(crs(test_rast)) %>%
#   terra::vect() %>%
#   terra::rasterize(y = test_rast,
#                    field = "pred_binary",
#                    fun = mean, na.rm = TRUE)
# 
# # get the extent of this particular raster, and crop it accordingly
# tempExt <- crds(plotObs_raw, na.rm = TRUE)
# 
# plotObs_raw_2 <- plotObs_raw %>%
#   crop(ext(min(tempExt[,1]), max(tempExt[,1]),
#            min(tempExt[,2]), max(tempExt[,2]))
#        )
# 
# map_ecoRegion_current_raw <- ggplot() +
#   geom_spatraster(data = plotObs_raw_2, maxcell = 50000000) +
#   geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
#   geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
# labs(title = c("Forest classification model predictions using
#                     contemporary dayMet data (from 2016)"),
#      subtitle = c("Values are the average binary predictions within each 8km x 8 km grid cell")) +
#   scale_fill_gradient2(low = "brown",
#                        mid = "wheat" ,
#                        high = "darkgreen" ,
#                        midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
#                        labels = c("grassShrub", "", "", "", "forest")) +
#   xlim(st_bbox(plotObs_2)[c(1,3)]) +
#   ylim(st_bbox(plotObs_2)[c(2,4)])

# for future v1
plotObs_future1 <- preds_future1_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "pred_binary", 
                   fun = mean, na.rm = TRUE)  
    
  
# get the extent of this particular raster, and crop it accordingly

plotObs_future1_2 <- plotObs_future1 %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_future1 <- ggplot() +
  geom_spatraster(data = plotObs_future1_2, maxcell = 50000000) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Forest classification predictions using future 
climate data from the BNU-ESM model"),
     subtitle = c("Values are the average binary predictions within each 8km x 8 km grid cell")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "", "", "", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])

# for future v2
plotObs_future2 <- preds_future2_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "pred_binary", 
                   fun = mean, na.rm = TRUE)  
    
# get the extent of this particular raster, and crop it accordingly
plotObs_future2_2 <- plotObs_future2 %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_future2 <- ggplot() +
  geom_spatraster(data = plotObs_future2_2, maxcell = 50000000) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future2_2)),fill=NA ) +
labs(title = paste0("Forest classification predictions using future 
climate data from the IPSL-CM5A-MR model")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future2_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future2_2)[c(2,4)])

## Plot the maps together
  ggarrange(map_ecoRegion_current, map_ecoRegion_future1, map_ecoRegion_future2, nrow = 1) %>% 
  annotate_figure(fig.lab = "Model Predictions of Forest Classification with Contemporary and Forecasted Climate Data", fig.lab.size = 20, 
                  bottom = text_grob("Note: Pink lines show the boundaries of contemporary grass/shrub and forest ecoregions, based on aggregated, EPA level 1 ecoregions"))

We used the binary predictions of forest occurrence shown above to identify the data we used to fit the ‘ecoregion-level’ cover models, but we also want to use the continuous predictions of forest probability to scale the cover model predictions to make them CONUS-wide. Because we use a “custom” threshold to convert the continuous probability predictions from the binomial model to binary classifications, we also need to scale the continuous predictions according to this custom threshold. We do this by performing linear interpolation of the continuous probability values above and below the threshold value, the results of which are shown below using both contemporary and forecasted climate data.

# setup linear interpolation functions
approxfun_pred <- approxfun(x = c(0, 0.3, 1), 
          y =c(0, 0.5, 1), 
          method = "linear")
# contemporary climate data
preds_byHand$preds_interp <- NA
preds_byHand[!is.na(preds_byHand$pred), "preds_interp"] <- approxfun_pred(preds_byHand[ !is.na(preds_byHand$pred), "pred"])

# predict with modeled future climate data v1
preds_future1_byHand$preds_interp <- NA
preds_future1_byHand[!is.na(preds_future1_byHand$pred), "preds_interp"] <- approxfun_pred(preds_future1_byHand[ !is.na(preds_future1_byHand$pred), "pred"])

# predict with modeled future climate data v2
preds_future2_byHand$preds_interp <- NA
preds_future2_byHand[!is.na(preds_future2_byHand$pred), "preds_interp"] <- approxfun_pred(preds_future2_byHand[ !is.na(preds_future2_byHand$pred), "pred"])

## Plot predictions with contemporary data
# make into a raster
plotObs <- preds_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "preds_interp", 
                   fun = mean, na.rm = TRUE) 

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)

plotObs_2 <- plotObs %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_current <- ggplot() +
  geom_spatraster(data = plotObs_2, maxcell = 50000000) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = c("Forest probability predictions (interpolated) using 
                    contemporary dayMet data (from 2016)"),
     subtitle = c("Values are the average binary predictions within each 8km x 8 km grid cell")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "", "", "", "forest")) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])


# for future v1
plotObs_future1 <- preds_future1_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "preds_interp", 
                   fun = mean, na.rm = TRUE)  
    
  
# get the extent of this particular raster, and crop it accordingly

plotObs_future1_2 <- plotObs_future1 %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_future1 <- ggplot() +
  geom_spatraster(data = plotObs_future1_2, maxcell = 50000000) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Forest probability predictions (interpolated) using future 
climate data from the BNU-ESM model"),
     subtitle = c("Values are the average binary predictions within each 8km x 8 km grid cell")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "", "", "", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])

# for future v2
plotObs_future2 <- preds_future2_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "preds_interp", 
                   fun = mean, na.rm = TRUE)  
    
# get the extent of this particular raster, and crop it accordingly
plotObs_future2_2 <- plotObs_future2 %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_future2 <- ggplot() +
  geom_spatraster(data = plotObs_future2_2, maxcell = 50000000) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future2_2)),fill=NA ) +
labs(title = paste0("Forest probability predictions (interpolated) using future 
climate data from the IPSL-CM5A-MR model")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future2_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future2_2)[c(2,4)])

## Plot the maps together
  ggarrange(map_ecoRegion_current, map_ecoRegion_future1, map_ecoRegion_future2, nrow = 1) %>% 
  annotate_figure(fig.lab = "Model Predictions of Continuous Forest Probability - interpolated- with Contemporary and Forecasted Climate Data", fig.lab.size = 20, 
                  bottom = text_grob("Note: Pink lines show the boundaries of contemporary grass/shrub and forest ecoregions, based on aggregated, EPA level 1 ecoregions"))

Now, we generate CONUS-wide model predictions of absolute cover

First, we use our models to predict absolute cover at the ecoregion or CONUS scale, depending on the functional type, using both contemporary climate and weather data, and then with future climate and weather data downscaled from the two GCMs. For some functional types we were able to fit a model that worked well across all of CONUS, while for other types, models performed better when we fit them to a single ecoregion at a time. Below is a list of the models we use:

• total herbaceous cover in the grass/shrub ecoregion
• total herbaceous cover in the forest ecoregion
• total tree cover in the grass/shrub ecoregion
• total tree cover in the forest ecoregion
• total shrub cover in CONUS
• total bare ground cover in CONUS
• The proportion of total tree cover that is needle-leaved in the grass/shrub ecoregion
• The proportion of total tree cover that is needle-leaved in the forest ecoregion
• The proportion of total tree cover that is broad-leaved in the grass/shrub ecoregion
• The proportion of total tree cover that is broad-leaved in the forest ecoregion
• The proportion of total herbaceous cover that is forbs in CONUS
• The proportion of total herbaceous cover that is C3 grass in CONUS
• The proportion of total herbaceous cover that is C4 grass in CONUS

*Note: Complete equations for each model used here can be found in the document “03_CompareModels.html”

# Total herbaceous cover: 
# Grass/shrub: 1SE lambda model
totalHerb_GS <-  readRDS("./models/betaLASSO/TotalHerbaceousCover_noTrees_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# Forest: best lambda model
totalHerb_F <-  readRDS("./models/betaLASSO/TotalHerbaceousCover_trees_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# Total tree cover: 
# Grass/shrub: best lambda model 
# totalTree_GS <- # include anomalies option
#   #readRDS("./models/betaLASSO/TotalTreeCover_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
#   #remove anomalies option #
#   readRDS("./models/betaLASSO/TotalTreeCover_trees_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest: best lambda model 
totalTree_F <- readRDS("./models/betaLASSO/TotalTreeCover_trees_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
  
# Bare ground: 
# CONUS - 1SE lambda model
bareGround_CONUS <-  readRDS("./models/betaLASSO/BareGroundCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")

# Shrub: 
# CONUS - best lambda model
shrub_CONUS <- readRDS("./models/betaLASSO/ShrubCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# Broad-leaved tree: 
# Grass/shrub: best lambda model
# BLtree_GS <- readRDS("./models/betaLASSO/AngioTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest:  best lambda model
BLtree_F <- readRDS("./models/betaLASSO/AngioTreeCover_prop_trees_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_halfSELambdaGLM.rds")
  
# Needle-leaved tree: 
# Grass/shrub: best lambda model
# NLtree_GS <-  readRDS("./models/betaLASSO/ConifTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest:  best lambda model
NLtree_F <- readRDS("./models/betaLASSO/ConifTreeCover_prop_trees_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# Forb:
# CONUS: best lambda model
forb_CONUS <-  readRDS("./models/betaLASSO/ForbCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# C3 grass:
# CONUS: 1SE lambda model
C3grass_CONUS <- readRDS("./models/betaLASSO/C3GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# C4 grass: 
# CONUS: best lambda model
C4grass_CONUS <- readRDS("./models/betaLASSO/C4GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

## Now, predict absolute cover with scaled climate/weather/soils data 
## with contemporary climate data

totalHerb_F_predsContemp <- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsContemp <- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsContemp <- climDatPred %>% 
  mutate(absTotalTree_GS = 0)
  # makePredictions(predictionDF = climDatPred,
  #                                           modelObject = totalTree_GS) %>% 
  # rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")

BLtree_GS_predsContemp<- climDatPred %>% 
  mutate(percBLTree_GS = 0)
  # makePredictions(predictionDF = climDatPred,
  #                                           modelObject = BLtree_GS) %>% 
  # rename(percBLTree_GS = "modelPreds")

NLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsContemp <- climDatPred %>% 
  mutate(percNLTree_GS = 0)
  # makePredictions(predictionDF = climDatPred,
  #                                           modelObject = NLtree_GS) %>% 
  # rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

contempPreds <- totalHerb_F_predsContemp  %>% 
  cbind(totalHerb_GS_predsContemp %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsContemp %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsContemp %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsContemp %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsContemp %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsContemp %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsContemp %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsContemp %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsContemp %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsContemp %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsContemp %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsContemp %>% select(percForb_CONUS)) %>% 
  cbind(preds_byHand %>% select(preds_interp)) %>% 
  rename(prob_Forest = "preds_interp") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

# get predictions with first climate model 
## with contemporary climate data
totalHerb_F_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture1 <- forecastClimSoilsDatPred_1 %>% 
  mutate(absTotalTree_GS = 0)
  # makePredictions(predictionDF = forecastClimSoilsDatPred_1,
  #                                           modelObject = totalTree_GS) %>% 
  # rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture1 <- forecastClimSoilsDatPred_1 %>% 
  mutate(percBLTree_GS = 0)
  # makePredictions(predictionDF = forecastClimSoilsDatPred_1,
  #                                           modelObject = BLtree_GS) %>% 
  # rename(percBLTree_GS = "modelPreds")

NLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture1 <- forecastClimSoilsDatPred_1 %>% 
  mutate(percNLTree_GS = 0)
  # makePredictions(predictionDF = forecastClimSoilsDatPred_1,
  #                                           modelObject = NLtree_GS) %>% 
  # rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

predsFuture1 <- totalHerb_F_predsFuture1  %>% 
  cbind(totalHerb_GS_predsFuture1 %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsFuture1 %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsFuture1 %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsFuture1 %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsFuture1 %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsFuture1 %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsFuture1 %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsFuture1 %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsFuture1 %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsFuture1 %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsFuture1 %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsFuture1 %>% select(percForb_CONUS)) %>% 
  cbind(preds_future1_byHand %>% select(preds_interp)) %>% 
  rename(prob_Forest = "preds_interp") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

# get predictions with second climate model
totalHerb_F_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture2 <- forecastClimSoilsDatPred_2 %>% 
  mutate(absTotalTree_GS = 0)
  # makePredictions(predictionDF = forecastClimSoilsDatPred_2,
  #                                           modelObject = totalTree_GS) %>% 
  # rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture2 <- forecastClimSoilsDatPred_2 %>% 
  mutate(percBLTree_GS = 0)
  # makePredictions(predictionDF = forecastClimSoilsDatPred_2,
  #                                           modelObject = BLtree_GS) %>% 
  # rename(percBLTree_GS = "modelPreds")

NLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture2 <- forecastClimSoilsDatPred_2 %>% 
  mutate(percNLTree_GS = 0)
  # makePredictions(predictionDF = forecastClimSoilsDatPred_2,
  #                                           modelObject = NLtree_GS) %>% 
  # rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

predsFuture2 <- totalHerb_F_predsFuture2  %>% 
  cbind(totalHerb_GS_predsFuture2 %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsFuture2 %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsFuture2 %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsFuture2 %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsFuture2 %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsFuture2 %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsFuture2 %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsFuture2 %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsFuture2 %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsFuture2 %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsFuture2 %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsFuture2 %>% select(percForb_CONUS)) %>% 
  cbind(preds_future2_byHand %>% select(preds_interp)) %>% 
  rename(prob_Forest = "preds_interp") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

Once we generate predictions using each of these models, we transform the ecoregion-level predictions for total herbaceous cover and total tree cover into CONUS-wide predictions by scaling them according to predicted ecoregion probability from the forest classification model. We retain the ecoregion-level predictions for later use in calculating CONUS-wide predictions of ‘level 2’ functional types

e.g.:absolute total tree cover in CONUS = (absolute total tree cover in grass/shrub * P(grass/shrub)) + (absolute total tree cover in forest * P(forest))

## for contemporary data
contempPreds <- contempPreds %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F)
         )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F)
         )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F)
         )

Then, for those models that are predicting proportions, we scale all of the relevant predictions so they sum to one. If models are ecoregion-level, we do this process at the ecoregion scale. (e.g. Model predictions of the proportion of total herbaceous cover that is forbs, the proportion of total herbaceous cover that is C3 grass, and the proportion of total herbaceous cover that is C4 grass should sum to 1, so: scaled proportion of total herbaceous cover that is forbs = proportion of total herbaceous cover that is forbs/(proportion of total herbaceous cover that is forbs + proportion of total herbaceous cover that is C3 grass + proportion of total herbaceous cover that is C4 grass))

contempPreds <- contempPreds %>% 
    mutate(percBLTree_scaled_GS = 0,#percBLTree_GS/(percBLTree_GS + percNLTree_GS),
         percNLTree_scaled_GS = 0,#percNLTree_GS/(percBLTree_GS + percNLTree_GS),
         percBLTree_scaled_F = percBLTree_F/(percBLTree_F + percNLTree_F),
         percNLTree_scaled_F = percNLTree_F/(percBLTree_F + percNLTree_F),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

# ggplot(contempPreds) +
#  geom_point(aes(x, y, col = percBLTree_GS)) #+
 # geom_density(aes(percBLTree_scaled_GS)) +
 # geom_density(aes(absTotalTree_F), col = "green")

predsFuture1 <- predsFuture1  %>% 
    mutate(percBLTree_scaled_GS = 0,#percBLTree_GS/(percBLTree_GS + percNLTree_GS),
         percNLTree_scaled_GS = 0,#percNLTree_GS/(percBLTree_GS + percNLTree_GS),
         percBLTree_scaled_F = percBLTree_F/(percBLTree_F + percNLTree_F),
         percNLTree_scaled_F = percNLTree_F/(percBLTree_F + percNLTree_F),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

predsFuture2 <- predsFuture2  %>% 
    mutate(percBLTree_scaled_GS = 0,#percBLTree_GS/(percBLTree_GS + percNLTree_GS),
         percNLTree_scaled_GS = 0,#percNLTree_GS/(percBLTree_GS + percNLTree_GS),
         percBLTree_scaled_F = percBLTree_F/(percBLTree_F + percNLTree_F),
         percNLTree_scaled_F = percNLTree_F/(percBLTree_F + percNLTree_F),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

Then, we convert these predicted proportions into predictions of absolute cover. For the functional types that are predicted at the ecoregion level, we do this process for each ecoregion. (e.g. absolute cover of needle-leaved trees in grass/shrub = scaled proportion of total tree cover that is needle-leaved in grass/shrub * absolute cover of total treesin grass/shrub)

## for contemporary data
contempPreds <- contempPreds %>% 
    mutate(absNLTree_GS = 0,#(percNLTree_scaled_GS * absTotalTree_GS), 
         absBLTree_GS = 0,#(percBLTree_scaled_GS * absTotalTree_GS), 
         absNLTree_F = (percNLTree_scaled_F * absTotalTree_F), 
         absBLTree_F = (percBLTree_scaled_F * absTotalTree_F), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C3 grass isn't ecoregion-level
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C4 grass isn't ecoregion-level
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS) # can use CONUS-wide total herbaceous since forb grass isn't ecoregion-level
  )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
    mutate(absNLTree_GS = 0,#(percNLTree_scaled_GS * absTotalTree_GS), 
         absBLTree_GS = 0,#(percBLTree_scaled_GS * absTotalTree_GS),
         absNLTree_F = (percNLTree_scaled_F * absTotalTree_F), 
         absBLTree_F = (percBLTree_scaled_F * absTotalTree_F), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C3 grass isn't ecoregion-level
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C4 grass isn't ecoregion-level
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS) # can use CONUS-wide total herbaceous since forb grass isn't ecoregion-level
  )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
    mutate(absNLTree_GS = 0,#(percNLTree_scaled_GS * absTotalTree_GS), 
         absBLTree_GS = 0,#(percBLTree_scaled_GS * absTotalTree_GS),
         absNLTree_F = (percNLTree_scaled_F * absTotalTree_F), 
         absBLTree_F = (percBLTree_scaled_F * absTotalTree_F), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C3 grass isn't ecoregion-level
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C4 grass isn't ecoregion-level
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS) # can use CONUS-wide total herbaceous since forb grass isn't ecoregion-level
  )

Finally, for those functional types that are predicted at the ecoregion level, we scale them to be CONUS-wide.

contempPreds <- contempPreds %>% 
  mutate(absBLTree_CONUS = (prob_grassShrub * absBLTree_GS) + (prob_Forest * absBLTree_F),
         absNLTree_CONUS = (prob_grassShrub * absNLTree_GS) + (prob_Forest * absNLTree_F)
  )

# ggplot(contempPreds) +
#  geom_point(aes(x, y, col = absNLTree_CONUS))

predsFuture1 <- predsFuture1 %>% 
    mutate(absBLTree_CONUS = (prob_grassShrub * absBLTree_GS) + (prob_Forest * absBLTree_F),
         absNLTree_CONUS = (prob_grassShrub * absNLTree_GS) + (prob_Forest * absNLTree_F)
  )

predsFuture2 <- predsFuture2 %>% 
    mutate(absBLTree_CONUS = (prob_grassShrub * absBLTree_GS) + (prob_Forest * absBLTree_F),
         absNLTree_CONUS = (prob_grassShrub * absNLTree_GS) + (prob_Forest * absNLTree_F)
  )

contempPreds <- contempPreds %>% 
  mutate(
     absNLTrees_V1 = (#(prob_grassShrub * percNLTree_GS) + 
       (prob_Forest * percNLTree_F))/(#(prob_grassShrub * percBLTree_GS) + 
         (prob_Forest * percBLTree_F) + (
           #(prob_grassShrub * percNLTree_GS) + 
           (prob_Forest * percNLTree_F))) * (#(prob_grassShrub * absTotalTree_GS) + 
             (prob_Forest * absTotalTree_F)) ,# option 1 (old approach)
         
         absNLTrees_V2 =  (prob_grassShrub * 0#((percNLTree_GS/(percBLTree_GS + percNLTree_GS)) * absTotalTree_GS)
                           ) + (prob_Forest * 
((percNLTree_F/(percBLTree_F + percNLTree_F)) * absTotalTree_F)# option 2 (new approach)
)
)
# 
# ggplot(contempPreds) + 
#   geom_density(aes(absNLTree_CONUS)) + # in-code calculation
#   geom_density(aes(absNLTrees_V1), col = "red") + # old version of calculation
#   geom_density(aes(absNLTrees_V2), col = "blue") # new version of calculation

Below are maps of predicted absolute cover for each functional group we modeled, after predictions have been scaled using the process outlined above. These maps show model-predicted cover using contemporary climate and weather data, residuals when comparing these contemporary predictions to observations, and changes in model predictions of cover when predictions are made using simulated future climate and weather data from two GCMs.

# prep observed data for use in figures
obsDat <- modDat_1_s %>% 
  select(x, y, Year, ShrubCover, TotalHerbaceousCover, TotalTreeCover, BareGroundCover, C3GramCover_prop, C4GramCover_prop, ForbCover_prop, AngioTreeCover_prop, ConifTreeCover_prop) %>% 
  mutate(absNLTree_CONUS = ConifTreeCover_prop * TotalTreeCover,
         absBLTree_CONUS = AngioTreeCover_prop * TotalTreeCover, 
         absC3grass_CONUS = C3GramCover_prop * TotalHerbaceousCover, 
         absC4grass_CONUS = C4GramCover_prop * TotalHerbaceousCover,
         absForb_CONUS = ForbCover_prop * TotalHerbaceousCover
         ) %>% 
  rename(absShrub_CONUS = ShrubCover,
         absTotalTree_CONUS = TotalTreeCover, 
         absTotalHerb_CONUS = TotalHerbaceousCover, 
         absBareGround_CONUS = BareGroundCover)

# list of absolute cover names
responseNames <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")

absoluteCoverMaps_contemp <- lapply(responseNames, 
                                      FUN = function(x) {
                                       ## contemporary absolute cover
                                        # rasterize response 
                                        resp_rast <- contempPreds %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_2 <- resp_rast %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # make a map
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of absolute cover of ",x,", 
                                                    \n after scaling according to ecoregion prediction")) +
                                          scale_fill_gradient2(low = "brown",
                                                               mid = "wheat" ,
                                                               high = "darkgreen" , 
                                                               midpoint = 0, limits = c(0,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of absolute cover of ",x)) +
                                          theme_classic()
                                        
                                        ## calculate residuals
                                        # make maps of residuals between absolute cover predictions and observations
                                        # rasterize observations 
                                        plotObservations <- obsDat %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs("EPSG:4326")) %>% 
                                          terra::project(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE) 
                                         plotObservations <- plotObservations %>%  
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate residuals (observed - predicted)
                                        resids_rast <- plotObservations - resp_rast_2
                                        ## aggregate the residuals raster to make it easier to see 
                                        resids_rastCoarse <- resids_rast %>% 
                                          aggregate(fact = 2, fun = "mean", na.rm = TRUE)
                                        # map the residuals
                                        
                                        mapResids <- ggplot() +
                                          geom_spatraster(data = resids_rastCoarse) + 
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% 
                                                    st_crop(st_bbox(resids_rast)),fill=NA ) +
                                          labs(title = paste0("Residuals (observations - predicted absolute cover) for ",x)) +
                                          scale_fill_gradient2(low = "red",
                                                               mid = "white" ,
                                                               high = "blue" , 
                                                               midpoint = 0,   na.value = "lightgrey",
                                                               limits = c(-1,1)
                                          )  + 
                                          xlim(st_bbox(resids_rast)[c(1,3)]) + 
                                          ylim(st_bbox(resids_rast)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram of residuals
                                        histResids <- ggplot() + 
                                          geom_density(aes(terra::values(resids_rast$mean, na.rm = TRUE, mat = FALSE)), 
                                                       col = "black", fill = "darkgrey") + 
                                          geom_vline(aes(xintercept = 0), linetype = 2) +
                                          xlab(paste0("Residuals (obs - pred) ",x)) +
                                          theme_classic()
                                        
                                         ## future model 1 absolute cover
                                        # rasterize response 
                                        resp_rast_future1 <- predsFuture1 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future1_2 <- resp_rast_future1 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
                                        # make a map
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of absolute cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        ## future model 2 absolute cover
                                        # rasterize response 
                                        resp_rast_future2 <- predsFuture2 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future2_2 <- resp_rast_future2 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
                                        # make a map
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of absolute cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        
                                        
                                        return(list("rast" = resp_rast_2, 
                                                    "map" = resp_map, 
                                                    "hist" = resp_hist,
                                                    "rast_resids" = resids_rast,
                                                    "map_resids" = mapResids,
                                                    "hist_resids" = histResids,
                                                    "rast_future1" = resp_rast_future1_2,
                                                    "map_future1" = resp_map_future1, 
                                                    "hist_future1" = resp_hist_future1,
                                                    "rast_future2" = resp_rast_future2_2,
                                                    "map_future2" = resp_map_future2, 
                                                    "hist_future2" = resp_hist_future2))
                                      })


names(absoluteCoverMaps_contemp) <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")

ggarrange(
ggarrange(absoluteCoverMaps_contemp$absTotalTree_CONUS$map, 
          absoluteCoverMaps_contemp$absTotalTree_CONUS$map_resids,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future1,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future2,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$hist,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_resids,
           absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absTotalHerb_CONUS$map, 
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_resids,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future1,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future2,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absShrub_CONUS$map, 
          absoluteCoverMaps_contemp$absShrub_CONUS$map_resids,
          absoluteCoverMaps_contemp$absShrub_CONUS$map_future1,
          absoluteCoverMaps_contemp$absShrub_CONUS$map_future2,
          absoluteCoverMaps_contemp$absShrub_CONUS$hist,
          absoluteCoverMaps_contemp$absShrub_CONUS$hist_resids,
          absoluteCoverMaps_contemp$absShrub_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absShrub_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absBareGround_CONUS$map, 
          absoluteCoverMaps_contemp$absBareGround_CONUS$map_resids,
          absoluteCoverMaps_contemp$absBareGround_CONUS$map_future1,
          absoluteCoverMaps_contemp$absBareGround_CONUS$map_future2,
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist,
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absNLTree_CONUS$map, 
          absoluteCoverMaps_contemp$absNLTree_CONUS$map_resids,
          absoluteCoverMaps_contemp$absNLTree_CONUS$map_future1,
          absoluteCoverMaps_contemp$absNLTree_CONUS$map_future2,
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist,
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absBLTree_CONUS$map, 
          absoluteCoverMaps_contemp$absBLTree_CONUS$map_resids,
          absoluteCoverMaps_contemp$absBLTree_CONUS$map_future1,
          absoluteCoverMaps_contemp$absBLTree_CONUS$map_future2,
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist,
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist_resids,
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absC3grass_CONUS$map, 
          absoluteCoverMaps_contemp$absC3grass_CONUS$map_resids,
          absoluteCoverMaps_contemp$absC3grass_CONUS$map_future1,
          absoluteCoverMaps_contemp$absC3grass_CONUS$map_future2,
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist,
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absC4grass_CONUS$map, 
          absoluteCoverMaps_contemp$absC4grass_CONUS$map_resids,
          absoluteCoverMaps_contemp$absC4grass_CONUS$map_future1,
          absoluteCoverMaps_contemp$absC4grass_CONUS$map_future2,
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist,
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absForb_CONUS$map, 
          absoluteCoverMaps_contemp$absForb_CONUS$map_resids,
          absoluteCoverMaps_contemp$absForb_CONUS$map_future1,
          absoluteCoverMaps_contemp$absForb_CONUS$map_future2,
          absoluteCoverMaps_contemp$absForb_CONUS$hist,
          absoluteCoverMaps_contemp$absForb_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absForb_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absForb_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ), 
ncol = 1
)

Now, we show quantile plots for model model predictions of absolute cover made with contemporary data. These plots compare final, scaled predictions of absolute cover to observed absolute cover. However, instead of showing raw data, which can be uninformative with large datasets such as ours, we show the average model predicted value or observed value within a given quantile of a given predictor variable. In this case, we are showing values within each percentile. We also show the inter-quantile range (IQR) of model-predicted and response values in each percentile, which appear as bands behind the points showing the mean values. Note that the predictor variables included in each quantile plot are the predictors that appear in all of the models used to generate a given absolute cover prediction (e.g. the predictors variables included in the quantile plot for absolute needle-leaved tree cover are the predictors that are included in any o fthe following models: total tree cover-grass/shrub, total tree cover-forest, prop. needle-leaved tree cover-grass/shrub, and prop. needle-leaved cover-forest). Note also that the predictor variables in these figures are scaled.

# predict forest classification
modDat_ecoregionFitQuantile <- modDat_1_s #%>% 
  # #rename("Long" = x, "Lat" = y) %>% 
  # dplyr::select(c(newRegion, t_warm, t_cold, prcp_wet, annWatDef, prcpTempCorr, isothermality, soilDepth, sand, coarse, carbon,
  #                                        x, y)) %>% 
  # mutate(newRegion = as.factor(newRegion)) 

# predict with contemporary climate data 
preds_byHand_quantile <- modDat_ecoregionFitQuantile %>% 
  mutate(pred =  as.numeric(1/ (1 + exp(-(-0.67478165 +  
                                            1.71345065*((prcp - 483.34337211) / 321.52170566) + 
                                            -0.01633932*((prcp_seasonality - 0.97338094) / 0.23365057) + 
                                            -0.91285243*((prcpTempCorr - 0.01681097) / 0.39882063) +  
                                            0.53377439*((annWetDegDays - 1808.61043539) / 1083.25298131) + 
                                            -0.60523936*((annWatDef_95 - 164.68611646) / 104.14892564) +  
                                            0.20578879*((coarse - 9.66120017) / 10.16883531) + 
                                            -0.61568228*((AWHC - 14.61190409) / 5.65765542) + 
                                            -0.31818631*I(((tmean - 11.0291275) / 4.89048541)^2) + 
                                            0.20104126*I(((prcp_seasonality - 0.97338094) / 0.23365057)^2) + 
                                            -0.89055646*I(((prcpTempCorr - 0.01681097) / 0.39882063)^2) +  
                                            0.14226917*I(((isothermality - 37.76731082) / 5.2256935)^2) + 
                                            -0.11810211*I(((clay - 19.59171509) / 8.51944876)^2) +  
                                            0.11278028*I(((AWHC - 14.61190409) / 5.65765542)^2) +  
                                            0.91656456*((annWetDegDays - 1808.61043539) / 1083.25298131)*((annWatDef_95 - 164.68611646) / 104.14892564) + 
                                            -0.29650844*((prcp_dry - 3.18622575) / 5.36400552)*((isothermality - 37.76731082) / 5.2256935) + 
                                            -0.02061048*((prcpTempCorr - 0.01681097) / 0.39882063)*((isothermality - 37.76731082) / 5.2256935) +
                                            -0.28504255*((isothermality - 37.76731082) / 5.2256935)*((tmean - 11.0291275) / 4.89048541) + 
                                            0.50549093*((prcpTempCorr - 0.01681097) / 0.39882063)*((prcp_dry - 3.18622575) / 5.36400552) +
                                            -0.25285633*((prcpTempCorr - 0.01681097) / 0.39882063)*((tmean - 11.0291275) / 4.89048541) + 
                                            -0.16290284*((AWHC - 14.61190409) / 5.65765542)*((clay - 19.59171509) / 8.51944876) +  
                                            0.09163249*((coarse - 9.66120017) / 10.16883531)*((clay - 19.59171509) / 8.51944876))))))
# convert the continuous probability predictions to binary using previously-identified threshold 
preds_byHand_quantile <- preds_byHand_quantile %>% 
  mutate(pred_binary = pred,
    pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
         pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))

# interpolate the continuous predictions
preds_byHand_quantile$preds_interp <- NA
preds_byHand_quantile[!is.na(preds_byHand_quantile$pred), "preds_interp"] <- approxfun_pred(preds_byHand_quantile[ !is.na(preds_byHand_quantile$pred), "pred"])

## now, predict for each functional type 
# rename the data.frame so that the scaled predictors have the correct names for the models to predict with 
scaledPrednames <- modDat_1_s %>% 
  select(tmin_s:AWHC_s) %>% 
  names()
modDat_quantile <- modDat_1_s %>% 
  select(ShrubCover, TotalHerbaceousCover, TotalTreeCover, C3GramCover_prop, C4GramCover_prop, 
         ForbCover_prop, AngioTreeCover_prop, ConifTreeCover_prop, BareGroundCover, x, y, tmin_s:AWHC_s) %>% 
  rename_with(~str_remove(string = .x, pattern = "_s$"), any_of(scaledPrednames))
    ## with the observed data, convert proportions to absolute cover
modDat_quantile <- modDat_quantile %>% 
  mutate(C3GramCover_abs = C3GramCover_prop * TotalHerbaceousCover, 
         C4GramCover_abs = C4GramCover_prop * TotalHerbaceousCover, 
         ForbCover_abs = ForbCover_prop * TotalHerbaceousCover, 
         NLTreeCover_abs = ConifTreeCover_prop * TotalTreeCover, 
         BLTreeCover_abs = AngioTreeCover_prop * TotalTreeCover)
## with contemporary climate data
totalHerb_F_predsQuantile <- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = totalHerb_F, scale = FALSE) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsQuantile <- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = totalHerb_GS, scale = FALSE) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = totalTree_F, scale = FALSE) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsQuantile <- modDat_quantile %>% 
  mutate(absTotalTree_GS = 0)#makePredictions(predictionDF = modDat_quantile,
                                            #modelObject = totalTree_GS, scale = FALSE) %>% 
  #rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = bareGround_CONUS, scale = FALSE) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = shrub_CONUS, scale = FALSE) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = BLtree_F, scale = FALSE) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsQuantile <- modDat_quantile %>% 
  mutate(percBLTree_GS = 0)#makePredictions(predictionDF = modDat_quantile,
                                            #modelObject = BLtree_GS, scale = FALSE) %>% 
 # rename(percBLTree_GS = "modelPreds")

NLtree_F_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = NLtree_F, scale = FALSE) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsQuantile <- modDat_quantile %>% 
  mutate(percNLTree_GS = 0)
  #makePredictions(predictionDF = modDat_quantile,
  #                                           modelObject = NLtree_GS, scale = FALSE) %>% 
  # rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = C3grass_CONUS, scale = FALSE) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = C4grass_CONUS, scale = FALSE) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
                                            modelObject = forb_CONUS, scale = FALSE) %>% 
  rename(percForb_CONUS = "modelPreds")

preds_quantile <- totalHerb_F_predsQuantile  %>% 
  cbind(totalHerb_GS_predsQuantile %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsQuantile %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsQuantile %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsQuantile %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsQuantile %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsQuantile %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsQuantile %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsQuantile %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsQuantile %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsQuantile %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsQuantile %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsQuantile %>% select(percForb_CONUS)) %>% 
  cbind(preds_byHand_quantile %>% select(preds_interp)) %>% 
  rename(prob_Forest = "preds_interp") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest) %>% 
  cbind(preds_byHand_quantile %>% select(newRegion))

# for total herbaceous and total tree, scale according to ecoregion
preds_quantile <- preds_quantile %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
         )

## For percentages (NL/BL tree and C4/C3/Forb), scale so that they sum to one (for C4/C3/Forb, already have CONUS-wide models) (for NL/BL tree, have versions for both grass/shrub and forest)
preds_quantile <- preds_quantile %>% 
  mutate(percBLTree_scaled_GS = 0,#percBLTree_GS/(percBLTree_GS + percNLTree_GS),
         percNLTree_scaled_GS = 0,#percNLTree_GS/(percBLTree_GS + percNLTree_GS),
         percBLTree_scaled_F = percBLTree_F/(percBLTree_F + percNLTree_F),
         percNLTree_scaled_F = percNLTree_F/(percBLTree_F + percNLTree_F),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

## Convert percentages for level 2 cover into absolute cover values
preds_quantile <- preds_quantile %>% 
  mutate(absNLTree_GS = (percNLTree_scaled_GS * absTotalTree_GS), 
         absBLTree_GS = (percBLTree_scaled_GS * absTotalTree_GS), 
         absNLTree_F = (percNLTree_scaled_F * absTotalTree_F), 
         absBLTree_F = (percBLTree_scaled_F * absTotalTree_F), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C3 grass isn't ecoregion-level
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS), # can use CONUS-wide total herbaceous since C4 grass isn't ecoregion-level
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS) # can use CONUS-wide total herbaceous since forb grass isn't ecoregion-level
  )

## For tree level 2 cover, scale by ecoregion probability prediction to be CONUS_wide
preds_quantile <- preds_quantile %>%
  mutate(absBLTree_CONUS = (prob_grassShrub * absBLTree_GS) + (prob_Forest * absBLTree_F),
         absNLTree_CONUS = (prob_grassShrub * absNLTree_GS) + (prob_Forest * absNLTree_F)
         )

## with the observed data, convert proportions to absolute cover
preds_quantile <- preds_quantile %>% 
  mutate(C3GramCover_abs = C3GramCover_prop * TotalHerbaceousCover, 
         C4GramCover_abs = C4GramCover_prop * TotalHerbaceousCover, 
         ForbCover_abs = ForbCover_prop * TotalHerbaceousCover, 
         NLTreeCover_abs = ConifTreeCover_prop * TotalTreeCover, 
         BLTreeCover_abs = AngioTreeCover_prop * TotalTreeCover)

# total trees
# predictions for absolute needle leaved tree cover
deciles_totTree <- predvars2deciles(df = preds_quantile %>% select(TotalTreeCover, absTotalTree_CONUS, "tmean", "prcp",  "prcp_dry", "isothermality", "AWHC",  "prcpTempCorr", "clay", "carbon" ,"coarse", "sand", "prcp_seasonality") %>% drop_na()
                                   , response_vars = c("TotalTreeCover", "absTotalTree_CONUS"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp",  "prcp_dry", "isothermality", "AWHC",  "prcpTempCorr", "clay", "carbon" ,"coarse", "sand", "prcp_seasonality"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

# NL trees
# predictions for absolute needle leaved tree cover
deciles_NLtree <- predvars2deciles(df = preds_quantile %>% select(NLTreeCover_abs, absNLTree_CONUS, tmean,  prcp,  prcp_dry , isothermality, AWHC, prcpTempCorr, clay,  carbon, coarse, sand ,  prcp_seasonality,   prcpTempCorr_anom , isothermality_anom) %>% drop_na()
                                   , response_vars = c("NLTreeCover_abs", "absNLTree_CONUS"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean",  "prcp",  "prcp_dry" , "isothermality", "AWHC", "prcpTempCorr", "clay",  "carbon", "coarse", "sand" ,  "prcp_seasonality",   "prcpTempCorr_anom" , "isothermality_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
# # by hand just to test
# decilesPlot_byHand_NLTree <- preds_quantile %>% 
#   select(NLTreeCover_abs, absNLTree_CONUS, tmean,  prcp,  prcp_dry , isothermality, AWHC, prcpTempCorr, clay,  carbon, coarse, sand ,  prcp_seasonality,   prcpTempCorr_anom , isothermality_anom) %>% 
#   drop_na() %>% 
#   pivot_longer(cols = c(tmean, prcp,  prcp_dry , isothermality, AWHC, prcpTempCorr, clay,  carbon, coarse, sand ,  prcp_seasonality,   prcpTempCorr_anom , isothermality_anom)) %>% 
#   ggplot() + 
#   facet_wrap(~name) + 
#   geom_point(aes(x = value, y = NLTreeCover_abs), alpha = .1, col = "grey") + # observed
#   geom_point(aes(x = value, y = absNLTree_CONUS), alpha = .1, col = "lightblue") + # predicted
#   geom_smooth(aes(x = value, y = NLTreeCover_abs), col = "black") + 
#   geom_smooth(aes(x = value, y = absNLTree_CONUS), col = "blue")
# 
#  preds_quantile %>% 
#   select(NLTreeCover_abs, absNLTree_CONUS, tmean,  prcp,  prcp_dry , isothermality, AWHC, prcpTempCorr, clay,  carbon, coarse, sand ,  prcp_seasonality,   prcpTempCorr_anom , isothermality_anom) %>% 
#   drop_na() %>% 
#    ggplot() +
#    geom_point(aes(x = absNLTree_CONUS, y = NLTreeCover_abs)) +
#    geom_abline(aes(slope = 1, intercept = 0), col = "red") +
#    geom_smooth(aes(x = absNLTree_CONUS, y = NLTreeCover_abs))
#   

# # BL trees
# unique(c(names(coefficients(totalTree_F)), names(coefficients(totalTree_GS)), names(coefficients(BLtree_GS)), names(coefficients(BLtree_F))))
deciles_BLtree <- predvars2deciles(df = preds_quantile %>% select("BLTreeCover_abs", "absBLTree_CONUS",  "tmean",  "prcp",  "prcpTempCorr", "sand", "coarse",  "AWHC"  , "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality",     
"prcpTempCorr_anom" , "prcp_dry" , "carbon" ) %>% drop_na(), 
response_vars = c("BLTreeCover_abs", "absBLTree_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c("tmean",  "prcp",  "prcpTempCorr", "sand", "coarse",  "AWHC"  , "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality",     
"prcpTempCorr_anom" , "prcp_dry" , "carbon" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
                         
# shrubs
# unique(c(names(coefficients(shrub_CONUS))))
deciles_shrub <- predvars2deciles(df = preds_quantile %>% select("ShrubCover", "absShrub_CONUS","prcp",  "prcp_seasonality",  "prcpTempCorr", "sand" ,"coarse",  "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean"  ) %>% drop_na(), 
                                  response_vars = c("ShrubCover", "absShrub_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c( "prcp",  "prcp_seasonality",  "prcpTempCorr", "sand" ,"coarse",  "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean" ), # for predictors, combine list of predictors for all models 
cut_points = seq(0, 1, 0.01)
                                     )

# bare ground
# unique(c(names(coefficients())))
deciles_bareGround <- predvars2deciles(df = preds_quantile %>% select("BareGroundCover", "absBareGround_CONUS", "tmean"  ,"prcpTempCorr"   , "isothermality"   , "annWetDegDays"  ,
"sand"   ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ) %>% drop_na(), 
                                       response_vars = c("BareGroundCover", "absBareGround_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c( "tmean"  ,"prcpTempCorr"   , "isothermality"   , "annWetDegDays"  ,
"sand"   ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
# total herbaceous cover
#unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS))))
deciles_totalHerb <- predvars2deciles(df = preds_quantile %>% select("TotalHerbaceousCover", "absTotalHerb_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom",   "tmean",  
                "prcp_anom", "clay", "carbon") %>% drop_na(),
                response_vars = c("TotalHerbaceousCover", "absTotalHerb_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom",   "tmean",  
                "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

# forbs
# unique(c(names(coefficients(totalTree_F)), names(coefficients(totalTree_GS)), names(coefficients(BLtree_GS)), names(coefficients(BLtree_F))))
deciles_forb <- predvars2deciles(df = preds_quantile %>% select("ForbCover_abs", "absForb_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay", 
                "carbon", "prcp_seasonality",  "annWetDegDays" ,  "prcp_seasonality_anom", "prcpTempCorr_anom" ,  "annWetDegDays_anom") %>% drop_na(), 
                                 response_vars = c("ForbCover_abs", "absForb_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay", 
                "carbon", "prcp_seasonality",  "annWetDegDays" ,  "prcp_seasonality_anom", "prcpTempCorr_anom" ,  "annWetDegDays_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

# C3 grass
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C3grass <- predvars2deciles(df = preds_quantile %>% select("C3GramCover_abs", "absC3grass_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand",  "coarse",  "AWHC" ,  "isothermality_anom",  
                "tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom",  "prcpTempCorr_anom",  "annWetDegDays_anom", "prcp_seasonality" ) %>% drop_na(), 
                                    response_vars = c("C3GramCover_abs", "absC3grass_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand",  "coarse",  "AWHC" ,  "isothermality_anom",  
                "tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom",  "prcpTempCorr_anom",  "annWetDegDays_anom", "prcp_seasonality" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
# C4 grass 
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C4grass <- predvars2deciles(df = preds_quantile %>% select("C4GramCover_abs", "absC4grass_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC",  "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na(), 
                                    response_vars = c("C4GramCover_abs", "absC4grass_CONUS"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC",  "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

## save prediction data for comparison purposes 
saveRDS(preds_quantile, "../../../Data_processed/CoverData/IntermediateAnalysisFiles/ModelPredictionDataWithObservations.rds")
## "raw" proportion of trees that are needle leaved
## subset the predictions so that we are only making quantile plots using data points where we observations for both total tree and proportion of tree that is needle-leaved and broad-leaved
preds_herb <- preds_quantile %>% 
  select(TotalHerbaceousCover, absTotalHerb_F, absTotalHerb_GS,
         #C3GramCover_prop, C4GramCover_prop, ForbCover_prop, 
         tmean:AWHC, newRegion, prob_Forest, prob_grassShrub) %>%
  drop_na()

#### total herbaceous ####
# for the entire CONUS using the GS model
deciles_totalHerb_GS_entire <- predvars2deciles(df = preds_herb %>% select(TotalHerbaceousCover, absTotalHerb_GS,
                                                                               "prcp", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC") %>% drop_na()
                                   , response_vars = c("TotalHerbaceousCover", "absTotalHerb_GS"), # name of observations, followed by name of predictions
                                   pred_vars = c("prcp", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01))

quantPlot_totalHerb_GS_entire <- decile_dotplot_pq(df = deciles_totalHerb_GS_entire, response= c("TotalHerbaceousCover", "absTotalHerb_GS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Herbaceous Cover - Grass/shrub model in CONUS")

quantPlot_totalHerb_GS_entire <- add_dotplot_inset(quantPlot_totalHerb_GS_entire, df = deciles_totalHerb_GS_entire, response = c("TotalHerbaceousCover", "absTotalHerb_GS"), dfRaw = preds_herb %>% select(TotalHerbaceousCover, absTotalHerb_GS,  "prcp", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


# for only GS (EPA defined) using the GS model

deciles_totalHerb_GS_GS <- predvars2deciles(df = preds_herb %>% filter(newRegion == "dryShrubGrass") %>% select(TotalHerbaceousCover, absTotalHerb_GS,
                                                                               "prcp", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC") %>% drop_na()
                                   , response_vars = c("TotalHerbaceousCover", "absTotalHerb_GS"), # name of observations, followed by name of predictions
                                   pred_vars = c("prcp", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_totalHerb_GS_GS <- decile_dotplot_pq(df = deciles_totalHerb_GS_GS, response= c("TotalHerbaceousCover", "absTotalHerb_GS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Herbaceous Cover - Grass/shrub model in Grass/shrub")

quantPlot_totalHerb_GS_GS <- add_dotplot_inset(quantPlot_totalHerb_GS_GS, df = deciles_totalHerb_GS_GS, response = c("TotalHerbaceousCover", "absTotalHerb_GS"), dfRaw = preds_herb %>% filter(newRegion == "dryShrubGrass") %>% select(TotalHerbaceousCover, absTotalHerb_GS,
                                                                               "prcp", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC")  %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


# for the entire CONUS using the F model
deciles_totalHerb_F_entire <- predvars2deciles(df = preds_herb %>% select(TotalHerbaceousCover, absTotalHerb_F,
                                                                               "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na()
                                   , response_vars = c("TotalHerbaceousCover", "absTotalHerb_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_totalHerb_F_entire <- decile_dotplot_pq(df = deciles_totalHerb_F_entire, response= c("TotalHerbaceousCover", "absTotalHerb_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Herbaceous Cover - Forest model in CONUS")

quantPlot_totalHerb_F_entire <- add_dotplot_inset(quantPlot_totalHerb_F_entire, df = deciles_totalHerb_F_entire, response = c("TotalHerbaceousCover", "absTotalHerb_F"), dfRaw = preds_herb %>% select(TotalHerbaceousCover, absTotalHerb_F,  "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (EPA defined) using the F model
deciles_totalHerb_F_F <- predvars2deciles(df = preds_herb %>% filter(newRegion != "dryShrubGrass") %>% select(TotalHerbaceousCover, absTotalHerb_F,"prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na()
                                   , response_vars = c("TotalHerbaceousCover", "absTotalHerb_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_totalHerb_F_F <- decile_dotplot_pq(df = deciles_totalHerb_F_F, response= c("TotalHerbaceousCover", "absTotalHerb_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Herbaceous Cover - Forest model in Forest")

quantPlot_totalHerb_F_F <- add_dotplot_inset(quantPlot_totalHerb_F_F, df = deciles_totalHerb_F_F, response = c("TotalHerbaceousCover", "absTotalHerb_F"), dfRaw = preds_herb %>% filter(newRegion != "dryShrubGrass") %>% select(TotalHerbaceousCover, absTotalHerb_F,
                                                                              "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon")  %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


## "raw" proportion of trees that are needle leaved
## subset the predictions so that we are only making quantile plots using data points where we observations for both total tree and proportion of tree that is needle-leaved and broad-leaved
preds_trees <- preds_quantile %>% 
  select(TotalTreeCover, absTotalTree_F, absTotalTree_GS,
         ConifTreeCover_prop, percNLTree_GS, percNLTree_F, percNLTree_scaled_F, percNLTree_scaled_GS, NLTreeCover_abs, absNLTree_CONUS, absNLTree_GS, absNLTree_F,
         AngioTreeCover_prop, percBLTree_GS, percBLTree_F, percBLTree_scaled_F, percBLTree_scaled_GS, BLTreeCover_abs, absBLTree_CONUS,
         absBLTree_GS, absBLTree_F,
         tmean:AWHC, newRegion, prob_Forest, prob_grassShrub) %>%
  drop_na()

#### total trees ####
# # for the entire CONUS using the GS model
# deciles_totalTree_GS_entire <- predvars2deciles(df = preds_trees %>% select(TotalTreeCover, absTotalTree_GS,
#                                                                                prcp, prcp_seasonality, sand ,AWHC) %>% drop_na()
#                                    , response_vars = c("TotalTreeCover", "absTotalTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("prcp", "prcp_seasonality","sand" ,"AWHC"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01))
# 
# quantPlot_totalTree_GS_entire <- decile_dotplot_pq(df = deciles_totalTree_GS_entire, response= c("TotalTreeCover", "absTotalTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute Total Tree Cover - Grass/shrub model in CONUS")
# 
# quantPlot_totalTree_GS_entire <- add_dotplot_inset(quantPlot_totalTree_GS_entire, df = deciles_totalTree_GS_entire, response = c("TotalTreeCover", "absTotalTree_GS"), dfRaw = preds_trees %>% select(TotalTreeCover, absTotalTree_GS,  prcp, prcp_seasonality, sand ,AWHC) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# 
# # for only GS (EPA defined) using the GS model
# 
# deciles_totalTree_GS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(TotalTreeCover, absTotalTree_GS,
#                                                                                prcp, prcp_seasonality, sand ,AWHC) %>% drop_na()
#                                    , response_vars = c("TotalTreeCover", "absTotalTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("prcp", "prcp_seasonality","sand" ,"AWHC"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_totalTree_GS_GS <- decile_dotplot_pq(df = deciles_totalTree_GS_GS, response= c("TotalTreeCover", "absTotalTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute Total Tree Cover - Grass/shrub model in Grass/shrub")
# 
# quantPlot_totalTree_GS_GS <- add_dotplot_inset(quantPlot_totalTree_GS_GS, df = deciles_totalTree_GS_GS, response = c("TotalTreeCover", "absTotalTree_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(TotalTreeCover, absTotalTree_GS,
#                                                                                prcp, prcp_seasonality, sand ,AWHC)  %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # for only GS (model-defined using the GS model
# deciles_totalTree_GS_GSprob <- predvars2deciles(df = preds_trees %>% filter(prob_Forest<0.5) %>% select(TotalTreeCover, absTotalTree_GS,
#                                                                                prcp, prcp_seasonality, sand ,AWHC) %>% drop_na()
#                                    , response_vars = c("TotalTreeCover", "absTotalTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("prcp", "prcp_seasonality","sand" ,"AWHC"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_totalTree_GS_GSprob <- decile_dotplot_pq(df = deciles_totalTree_GS_GSprob, response= c("TotalTreeCover", "absTotalTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute Total Tree Cover - Grass/shrub model in Grass/shrub")
# 
# quantPlot_totalTree_GS_GSprob <- add_dotplot_inset(quantPlot_totalTree_GS_GSprob, df = deciles_totalTree_GS_GSprob, response = c("TotalTreeCover", "absTotalTree_GS"), dfRaw = preds_trees %>% filter(prob_Forest<0.5) %>% select(TotalTreeCover, absTotalTree_GS,
#                                                                                prcp, prcp_seasonality, sand ,AWHC)  %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for the entire CONUS using the F model
deciles_totalTree_F_entire <- predvars2deciles(df = preds_trees %>% select(TotalTreeCover, absTotalTree_F,
                                                                               tmean, prcp, prcp_dry, isothermality, AWHC, prcpTempCorr, clay, carbon, coarse, sand) %>% drop_na()
                                   , response_vars = c("TotalTreeCover", "absTotalTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp", "prcp_dry", "isothermality", "AWHC", "prcpTempCorr", "clay", "carbon", "coarse", "sand"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_totalTree_F_entire <- decile_dotplot_pq(df = deciles_totalTree_F_entire, response= c("TotalTreeCover", "absTotalTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Tree Cover - Forest model in CONUS")

quantPlot_totalTree_F_entire <- add_dotplot_inset(quantPlot_totalTree_F_entire, df = deciles_totalTree_F_entire, response = c("TotalTreeCover", "absTotalTree_F"), dfRaw = preds_trees %>% select(TotalTreeCover, absTotalTree_F,  tmean, prcp, prcp_dry, isothermality, AWHC, prcpTempCorr, clay, carbon, coarse, sand) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (EPA defined) using the F model
deciles_totalTree_F_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(TotalTreeCover, absTotalTree_F,tmean, prcp, prcp_dry, isothermality, AWHC, prcpTempCorr, clay, carbon, coarse, sand) %>% drop_na()
                                   , response_vars = c("TotalTreeCover", "absTotalTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp", "prcp_dry", "isothermality", "AWHC", "prcpTempCorr", "clay", "carbon", "coarse", "sand"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_totalTree_F_F <- decile_dotplot_pq(df = deciles_totalTree_F_F, response= c("TotalTreeCover", "absTotalTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Tree Cover - Forest model in Forest")

quantPlot_totalTree_F_F <- add_dotplot_inset(quantPlot_totalTree_F_F, df = deciles_totalTree_F_F, response = c("TotalTreeCover", "absTotalTree_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(TotalTreeCover, absTotalTree_F,
                                                                              tmean, prcp, prcp_dry, isothermality, AWHC, prcpTempCorr, clay, carbon, coarse, sand)  %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (model-defined using the F model
deciles_totalTree_F_Fprob <- predvars2deciles(df = preds_trees %>% filter(prob_Forest>=0.5) %>% select(TotalTreeCover, absTotalTree_F,
                                                                               tmean, prcp, prcp_dry, isothermality, AWHC, prcpTempCorr, clay, carbon, coarse, sand) %>% drop_na()
                                   , response_vars = c("TotalTreeCover", "absTotalTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp", "prcp_dry", "isothermality", "AWHC", "prcpTempCorr", "clay", "carbon", "coarse", "sand"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_totalTree_F_Fprob <- decile_dotplot_pq(df = deciles_totalTree_F_Fprob, response= c("TotalTreeCover", "absTotalTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Tree Cover - Forest model in Forest")

quantPlot_totalTree_F_Fprob <- add_dotplot_inset(quantPlot_totalTree_F_Fprob, df = deciles_totalTree_F_Fprob, response = c("TotalTreeCover", "absTotalTree_F"), dfRaw = preds_trees %>% filter(prob_Forest>=0.5) %>% select(TotalTreeCover, absTotalTree_F,
                                                                               prcp, prcp_seasonality, sand ,AWHC)  %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

#### NL trees ####
# # for the entire CONUS using the GS model
# deciles_NLtree_RawGS_entire <- predvars2deciles(df = preds_trees %>% select(ConifTreeCover_prop, percNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na()
#                                    , response_vars = c("ConifTreeCover_prop", "percNLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_RawGS_entire <- decile_dotplot_pq(df = deciles_NLtree_RawGS_entire, response= c("ConifTreeCover_prop", "percNLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Needle-Leaved - Grass/shrub model in CONUS")
# 
# quantPlot_NLtree_RawGS_entire <- add_dotplot_inset(quantPlot_NLtree_RawGS_entire, df = deciles_NLtree_RawGS_entire, response = c("ConifTreeCover_prop", "percNLTree_GS"), dfRaw = preds_trees %>% select(ConifTreeCover_prop, percNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # for only GS using the GS model (EPA-defined)
# deciles_NLtree_RawGS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na() 
#                                    , response_vars = c("ConifTreeCover_prop", "percNLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_RawGS_GS <- decile_dotplot_pq(df = deciles_NLtree_RawGS_GS, response= c("ConifTreeCover_prop", "percNLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Needle-Leaved - Grass/shrub model in Grass/shrub")
# 
# quantPlot_NLtree_RawGS_GS <- add_dotplot_inset(quantPlot_NLtree_RawGS_GS, df = deciles_NLtree_RawGS_GS, response = c("ConifTreeCover_prop", "percNLTree_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # for only GS using the GS model (model-defined)
# deciles_NLtree_RawGS_GSprob <- predvars2deciles(df = preds_trees %>% filter(prob_Forest < 0.5) %>% select(ConifTreeCover_prop, percNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na() 
#                                    , response_vars = c("ConifTreeCover_prop", "percNLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_RawGS_GSprob <- decile_dotplot_pq(df = deciles_NLtree_RawGS_GSprob, response= c("ConifTreeCover_prop", "percNLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Needle-Leaved - Grass/shrub model in Grass/shrub")
# 
# quantPlot_NLtree_RawGS_GSprob <- add_dotplot_inset(quantPlot_NLtree_RawGS_GSprob, df = deciles_NLtree_RawGS_GSprob, response = c("ConifTreeCover_prop", "percNLTree_GS"), dfRaw = preds_trees %>% filter(prob_Forest < 0.5) %>% select(ConifTreeCover_prop, percNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


# for the entire CONUS using the F model
deciles_NLtree_RawF_entire <- predvars2deciles(df = preds_trees %>% select(ConifTreeCover_prop, percNLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na()
                                   , response_vars = c("ConifTreeCover_prop", "percNLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_RawF_entire <- decile_dotplot_pq(df = deciles_NLtree_RawF_entire, response= c("ConifTreeCover_prop", "percNLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Needle-Leaved - forest model in CONUS")

quantPlot_NLtree_RawF_entire <- add_dotplot_inset(quantPlot_NLtree_RawF_entire, df = deciles_NLtree_RawF_entire, response = c("ConifTreeCover_prop", "percNLTree_F"), dfRaw = preds_trees %>% select(ConifTreeCover_prop, percNLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (EPA defined) using the F model
deciles_NLtree_RawF_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na() 
                                   , response_vars = c("ConifTreeCover_prop", "percNLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_RawF_F <- decile_dotplot_pq(df = deciles_NLtree_RawF_F, response= c("ConifTreeCover_prop", "percNLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Needle-Leaved - forest model in forest")

quantPlot_NLtree_RawF_F <- add_dotplot_inset(quantPlot_NLtree_RawF_F, df = deciles_NLtree_RawF_F, response = c("ConifTreeCover_prop", "percNLTree_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass")%>% select(ConifTreeCover_prop, percNLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (Model-defined) using the F model
deciles_NLtree_RawF_Fprob <- predvars2deciles(df = preds_trees %>% filter(prob_Forest > 0.5) %>% select(ConifTreeCover_prop, percNLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na() 
                                   , response_vars = c("ConifTreeCover_prop", "percNLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_RawF_Fprob <- decile_dotplot_pq(df = deciles_NLtree_RawF_Fprob, response= c("ConifTreeCover_prop", "percNLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Needle-Leaved - forest model in forest")

quantPlot_NLtree_RawF_Fprob <- add_dotplot_inset(quantPlot_NLtree_RawF_Fprob, df = deciles_NLtree_RawF_Fprob, response = c("ConifTreeCover_prop", "percNLTree_F"), dfRaw = preds_trees %>% filter(prob_Forest > 0.5) %>% select(ConifTreeCover_prop, percNLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

## prediction of proportion of NL tree in Grass/Shrub, scaled w/ BL tree to sum to 1 
# # grass/shrub model predicting in Grass/Shrub (EPA-defined)
# deciles_NLtree_scaledGS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_scaled_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na() 
#                                    , response_vars = c("ConifTreeCover_prop", "percNLTree_scaled_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_scaledGS_GS <- decile_dotplot_pq(df = deciles_NLtree_scaledGS_GS, response= c("ConifTreeCover_prop", "percNLTree_scaled_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Needle-Leaved \n(scaled to sum to 1 w/ BL tree) - Grass/shrub model in Grass/shrub")
# 
# quantPlot_NLtree_scaledGS_GS <- add_dotplot_inset(quantPlot_NLtree_scaledGS_GS, df = deciles_NLtree_scaledGS_GS, response = c("ConifTreeCover_prop", "percNLTree_scaled_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_scaled_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # grass/shrub model predicting in CONUS
# deciles_NLtree_scaledGS_entire <- predvars2deciles(df = preds_trees %>% select(ConifTreeCover_prop, percNLTree_scaled_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na()
#                                    , response_vars = c("ConifTreeCover_prop", "percNLTree_scaled_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_scaledGS_entire <- decile_dotplot_pq(df = deciles_NLtree_scaledGS_entire, response= c("ConifTreeCover_prop", "percNLTree_scaled_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Needle-Leaved \n(scaled to sum to 1 w/ BL tree) - Grass/shrub model in CONUS")
# 
# quantPlot_NLtree_scaledGS_entire <- add_dotplot_inset(quantPlot_NLtree_scaledGS_entire, df = deciles_NLtree_scaledGS_entire, response = c("ConifTreeCover_prop", "percNLTree_scaled_GS"), dfRaw = preds_trees %>% select(ConifTreeCover_prop, percNLTree_scaled_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
#  
# ## prediction of absolute cover of NL tree in Grass/Shrub
# # grass/shrub model predicting in Grass/Shrub (EPA-defined)
# deciles_NLtree_absoluteGS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(NLTreeCover_abs, absNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na() 
#                                    , response_vars = c("NLTreeCover_abs", "absNLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_absoluteGS_GS <- decile_dotplot_pq(df = deciles_NLtree_absoluteGS_GS, response= c("NLTreeCover_abs", "absNLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute cover of needle-leaved trees - \nGrass/shrub model in Grass/shrub")
# 
# quantPlot_NLtree_absoluteGS_GS <- add_dotplot_inset(quantPlot_NLtree_absoluteGS_GS, df = deciles_NLtree_absoluteGS_GS, response = c("NLTreeCover_abs", "absNLTree_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(NLTreeCover_abs, absNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # grass/shrub model predicting in CONUS
# deciles_NLtree_absoluteGS_entire <- predvars2deciles(df = preds_trees %>% select(NLTreeCover_abs, absNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na()
#                                    , response_vars = c("NLTreeCover_abs", "absNLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "AWHC", "prcpTempCorr", "isothermality", "prcpTempCorr_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_NLtree_absoluteGS_entire <- decile_dotplot_pq(df = deciles_NLtree_absoluteGS_entire, response= c("NLTreeCover_abs", "absNLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute cover of needle-leaved trees -\nGrass/shrub model in CONUS")
# 
# quantPlot_NLtree_absoluteGS_entire <- add_dotplot_inset(quantPlot_NLtree_absoluteGS_entire, df = deciles_NLtree_absoluteGS_entire, response = c("NLTreeCover_abs", "absNLTree_GS"), dfRaw = preds_trees %>% select(NLTreeCover_abs, absNLTree_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

## prediction of proportion of NL tree in forest, scaled w/ BL tree to sum to 1 
# forest model predicting in forest (EPA-defined)
deciles_NLtree_scaledF_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_scaled_F,"tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na() 
                                   , response_vars = c("ConifTreeCover_prop", "percNLTree_scaled_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_scaledF_F <- decile_dotplot_pq(df = deciles_NLtree_scaledF_F, response= c("ConifTreeCover_prop", "percNLTree_scaled_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Needle-Leaved \n(scaled to sum to 1 w/ BL tree) - forest model in forest")

quantPlot_NLtree_scaledF_F <- add_dotplot_inset(quantPlot_NLtree_scaledF_F, df = deciles_NLtree_scaledF_F, response = c("ConifTreeCover_prop", "percNLTree_scaled_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(ConifTreeCover_prop, percNLTree_scaled_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# forest model predicting in CONUS
deciles_NLtree_scaledF_entire <- predvars2deciles(df = preds_trees %>% select(ConifTreeCover_prop, percNLTree_scaled_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na()
                                   , response_vars = c("ConifTreeCover_prop", "percNLTree_scaled_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_scaledF_entire <- decile_dotplot_pq(df = deciles_NLtree_scaledF_entire, response= c("ConifTreeCover_prop", "percNLTree_scaled_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Needle-Leaved \n(scaled to sum to 1 w/ BL tree) - forest model in CONUS")

quantPlot_NLtree_scaledF_entire <- add_dotplot_inset(quantPlot_NLtree_scaledF_entire, df = deciles_NLtree_scaledF_entire, response = c("ConifTreeCover_prop", "percNLTree_scaled_F"), dfRaw = preds_trees %>% select(ConifTreeCover_prop, percNLTree_scaled_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

 
## prediction of absolute cover of NL tree in forest
# forest model predicting in forest (EPA-defined)
deciles_NLtree_absoluteF_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(NLTreeCover_abs, absNLTree_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na() 
                                   , response_vars = c("NLTreeCover_abs", "absNLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_absoluteF_F <- decile_dotplot_pq(df = deciles_NLtree_absoluteF_F, response= c("NLTreeCover_abs", "absNLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute cover of needle-leaved trees - \nforest model in forest")

quantPlot_NLtree_absoluteF_F <- add_dotplot_inset(quantPlot_NLtree_absoluteF_F, df = deciles_NLtree_absoluteF_F, response = c("NLTreeCover_abs", "absNLTree_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(NLTreeCover_abs, absNLTree_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# forest model predicting in CONUS
deciles_NLtree_absoluteF_entire <- predvars2deciles(df = preds_trees %>% select(NLTreeCover_abs, absNLTree_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na()
                                   , response_vars = c("NLTreeCover_abs", "absNLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_absoluteF_entire <- decile_dotplot_pq(df = deciles_NLtree_absoluteF_entire, response= c("NLTreeCover_abs", "absNLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute cover of needle-leaved trees - \nforest model in CONUS")

quantPlot_NLtree_absoluteF_entire <- add_dotplot_inset(quantPlot_NLtree_absoluteF_entire, df = deciles_NLtree_absoluteF_entire, response = c("NLTreeCover_abs", "absNLTree_F"), dfRaw = preds_trees %>% select(NLTreeCover_abs, absNLTree_F, "tmean", "prcp_dry", "sand", "carbon", "AWHC", "isothermality", "prcpTempCorr", "prcpTempCorr_anom", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


## prediction of absolute cover of NL tree, scaled by ecoregion probability to be CONUS-wide 
deciles_NLtree_absolute_CONUS <- predvars2deciles(df = preds_trees %>% select(NLTreeCover_abs, absNLTree_CONUS, AWHC, carbon, clay, coarse, isothermality_anom, isothermality, tmean, prcp, prcp_dry, prcp_seasonality, prcpTempCorr_anom, prcpTempCorr, sand) %>% drop_na()
                                   , response_vars = c("NLTreeCover_abs", "absNLTree_CONUS"), # name of observations, followed by name of predictions
                                   pred_vars = c("AWHC", "carbon", "clay", "coarse", "isothermality_anom", "isothermality", "tmean", "prcp", "prcp_dry", "prcp_seasonality", "prcpTempCorr_anom", "prcpTempCorr", "sand"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc NL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_NLtree_absolute_CONUS <- decile_dotplot_pq(df = deciles_NLtree_absolute_CONUS, response= c("NLTreeCover_abs", "absNLTree_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute cover of needle-leaved trees - CONUS-wide")

quantPlot_NLtree_absolute_CONUS <- add_dotplot_inset(quantPlot_NLtree_absolute_CONUS, df = deciles_NLtree_absolute_CONUS, response = c("NLTreeCover_abs", "absNLTree_CONUS"), dfRaw = preds_trees %>% select(NLTreeCover_abs, absNLTree_CONUS, "AWHC", "carbon", "clay", "coarse", "isothermality_anom", "isothermality", "tmean", "prcp", "prcp_dry", "prcp_seasonality", "prcpTempCorr_anom", "prcpTempCorr", "sand") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


#### BL trees ####
# for the entire CONUS using the GS model
# deciles_BLtree_RawGS_entire <- predvars2deciles(df = preds_trees %>% select(AngioTreeCover_prop, percBLTree_GS, tmean, prcp, prcpTempCorr, sand, coarse, AWHC, prcp_seasonality_anom, isothermality_anom, isothermality, prcp_seasonality, prcpTempCorr_anom) %>% drop_na()
#                                    , response_vars = c("AngioTreeCover_prop", "percBLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_RawGS_entire <- decile_dotplot_pq(df = deciles_BLtree_RawGS_entire, response= c("AngioTreeCover_prop", "percBLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Broad-leaved - Grass/shrub model in CONUS")
# 
# quantPlot_BLtree_RawGS_entire <- add_dotplot_inset(quantPlot_BLtree_RawGS_entire, df = preds_trees %>% select(AngioTreeCover_prop, percBLTree_GS, tmean, prcp, prcpTempCorr, sand, coarse, AWHC, prcp_seasonality_anom, isothermality_anom, isothermality, prcp_seasonality, prcpTempCorr_anom) %>% drop_na(), response = c("AngioTreeCover_prop", "percBLTree_GS"), dfRaw = preds_trees %>% select(AngioTreeCover_prop, percBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# 
# # for oBLy GS using the GS model (EPA-defined)
# deciles_BLtree_RawGS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na() 
#                                    , response_vars = c("AngioTreeCover_prop", "percBLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_RawGS_GS <- decile_dotplot_pq(df = deciles_BLtree_RawGS_GS, response= c("AngioTreeCover_prop", "percBLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Broad-leaved - Grass/shrub model in Grass/shrub")
# 
# quantPlot_BLtree_RawGS_GS <- add_dotplot_inset(quantPlot_BLtree_RawGS_GS, df = deciles_BLtree_RawGS_GS, response = c("AngioTreeCover_prop", "percBLTree_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # for oBLy GS using the GS model (model-defined)
# deciles_BLtree_RawGS_GSprob <- predvars2deciles(df = preds_trees %>% filter(prob_Forest < 0.5) %>% select(AngioTreeCover_prop, percBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na() 
#                                    , response_vars = c("AngioTreeCover_prop", "percBLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_RawGS_GSprob <- decile_dotplot_pq(df = deciles_BLtree_RawGS_GSprob, response= c("AngioTreeCover_prop", "percBLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Broad-leaved - Grass/shrub model in Grass/shrub")
# 
# quantPlot_BLtree_RawGS_GSprob <- add_dotplot_inset(quantPlot_BLtree_RawGS_GSprob, df = deciles_BLtree_RawGS_GSprob, response = c("AngioTreeCover_prop", "percBLTree_GS"), dfRaw = preds_trees %>% filter(prob_Forest < 0.5) %>% select(AngioTreeCover_prop, percBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 

# for the entire CONUS using the F model
deciles_BLtree_RawF_entire <- predvars2deciles(df = preds_trees %>% select(AngioTreeCover_prop, percBLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na()
                                   , response_vars = c("AngioTreeCover_prop", "percBLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), 
                                   # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_RawF_entire <- decile_dotplot_pq(df = deciles_BLtree_RawF_entire, response= c("AngioTreeCover_prop", "percBLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Broad-leaved - forest model in CONUS")

quantPlot_BLtree_RawF_entire <- add_dotplot_inset(quantPlot_BLtree_RawF_entire, df = deciles_BLtree_RawF_entire, response = c("AngioTreeCover_prop", "percBLTree_F"), dfRaw = preds_trees %>% select(AngioTreeCover_prop, percBLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (EPA defined) using the F model
deciles_BLtree_RawF_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na() 
                                   , response_vars = c("AngioTreeCover_prop", "percBLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_RawF_F <- decile_dotplot_pq(df = deciles_BLtree_RawF_F, response= c("AngioTreeCover_prop", "percBLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Broad-leaved - forest model in forest")

quantPlot_BLtree_RawF_F <- add_dotplot_inset(quantPlot_BLtree_RawF_F, df = deciles_BLtree_RawF_F, response = c("AngioTreeCover_prop", "percBLTree_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass")%>% select(AngioTreeCover_prop, percBLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# for only F (Model-defined) using the F model
deciles_BLtree_RawF_Fprob <- predvars2deciles(df = preds_trees %>% filter(prob_Forest > 0.5) %>% select(AngioTreeCover_prop, percBLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na() 
                                   , response_vars = c("AngioTreeCover_prop", "percBLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_RawF_Fprob <- decile_dotplot_pq(df = deciles_BLtree_RawF_Fprob, response= c("AngioTreeCover_prop", "percBLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Broad-leaved - forest model in forest")

quantPlot_BLtree_RawF_Fprob <- add_dotplot_inset(quantPlot_BLtree_RawF_Fprob, df = deciles_BLtree_RawF_Fprob, response = c("AngioTreeCover_prop", "percBLTree_F"), dfRaw = preds_trees %>% filter(prob_Forest > 0.5) %>% select(AngioTreeCover_prop, percBLTree_F, tmean, prcp_dry, sand, carbon, AWHC, isothermality, prcpTempCorr, prcpTempCorr_anom, isothermality_anom, coarse) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# ## prediction of proportion of BL tree in Grass/Shrub, scaled w/ BL tree to sum to 1 
# # grass/shrub model predicting in Grass/Shrub (EPA-defined)
# deciles_BLtree_scaledGS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_scaled_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na() 
#                                    , response_vars = c("AngioTreeCover_prop", "percBLTree_scaled_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_scaledGS_GS <- decile_dotplot_pq(df = deciles_BLtree_scaledGS_GS, response= c("AngioTreeCover_prop", "percBLTree_scaled_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Broad-leaved (scaled to sum to 1 w/ BL tree) - Grass/shrub model in Grass/shrub")
# 
# quantPlot_BLtree_scaledGS_GS <- add_dotplot_inset(quantPlot_BLtree_scaledGS_GS, df = deciles_BLtree_scaledGS_GS, response = c("AngioTreeCover_prop", "percBLTree_scaled_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_scaled_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # grass/shrub model predicting in CONUS
# deciles_BLtree_scaledGS_entire <- predvars2deciles(df = preds_trees %>% select(AngioTreeCover_prop, percBLTree_scaled_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na()
#                                    , response_vars = c("AngioTreeCover_prop", "percBLTree_scaled_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_scaledGS_entire <- decile_dotplot_pq(df = deciles_BLtree_scaledGS_entire, response= c("AngioTreeCover_prop", "percBLTree_scaled_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Proportion of trees that are Broad-leaved (scaled to sum to 1 w/ BL tree) - Grass/shrub model in CONUS")
# 
# quantPlot_BLtree_scaledGS_entire <- add_dotplot_inset(quantPlot_BLtree_scaledGS_entire, df = deciles_BLtree_scaledGS_entire, response = c("AngioTreeCover_prop", "percBLTree_scaled_GS"), dfRaw = preds_trees %>% select(AngioTreeCover_prop, percBLTree_scaled_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

 
## prediction of absolute cover of BL tree in Grass/Shrub
# grass/shrub model predicting in Grass/Shrub (EPA-defined)
# deciles_BLtree_absoluteGS_GS <- predvars2deciles(df = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(BLTreeCover_abs, absBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na() 
#                                    , response_vars = c("BLTreeCover_abs", "absBLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_absoluteGS_GS <- decile_dotplot_pq(df = deciles_BLtree_absoluteGS_GS, response= c("BLTreeCover_abs", "absBLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute cover of Broad-leaved trees - Grass/shrub model in Grass/shrub")
# 
# quantPlot_BLtree_absoluteGS_GS <- add_dotplot_inset(quantPlot_BLtree_absoluteGS_GS, df = deciles_BLtree_absoluteGS_GS, response = c("BLTreeCover_abs", "absBLTree_GS"), dfRaw = preds_trees %>% filter(newRegion == "dryShrubGrass") %>% select(BLTreeCover_abs, absBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)
# 
# # grass/shrub model predicting in CONUS
# deciles_BLtree_absoluteGS_entire <- predvars2deciles(df = preds_trees %>% select(BLTreeCover_abs, absBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na()
#                                    , response_vars = c("BLTreeCover_abs", "absBLTree_GS"), # name of observations, followed by name of predictions
#                                    pred_vars = c("tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
# cut_points = seq(0, 1, 0.01)
#                                      )
# quantPlot_BLtree_absoluteGS_entire <- decile_dotplot_pq(df = deciles_BLtree_absoluteGS_entire, response= c("BLTreeCover_abs", "absBLTree_GS"), IQR = TRUE,
#                         CI = FALSE
#                         ) + ggtitle("Absolute cover of Broad-leaved trees - Grass/shrub model in CONUS")
# 
# quantPlot_BLtree_absoluteGS_entire <- add_dotplot_inset(quantPlot_BLtree_absoluteGS_entire, df = deciles_BLtree_absoluteGS_entire, response = c("BLTreeCover_abs", "absBLTree_GS"), dfRaw = preds_trees %>% select(BLTreeCover_abs, absBLTree_GS, "tmean", "prcp", "prcpTempCorr", "sand", "coarse", "AWHC", "prcp_seasonality_anom", "isothermality_anom", "isothermality", "prcp_seasonality", "prcpTempCorr_anom") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

## prediction of proportion of BL tree in forest, scaled w/ BL tree to sum to 1 
# forest model predicting in forest (EPA-defined)
deciles_BLtree_scaledF_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_scaled_F,"tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na() 
                                   , response_vars = c("AngioTreeCover_prop", "percBLTree_scaled_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_scaledF_F <- decile_dotplot_pq(df = deciles_BLtree_scaledF_F, response= c("AngioTreeCover_prop", "percBLTree_scaled_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Broad-leaved (scaled to sum to 1 w/ BL tree) - forest model in forest")

quantPlot_BLtree_scaledF_F <- add_dotplot_inset(quantPlot_BLtree_scaledF_F, df = deciles_BLtree_scaledF_F, response = c("AngioTreeCover_prop", "percBLTree_scaled_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(AngioTreeCover_prop, percBLTree_scaled_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# forest model predicting in CONUS
deciles_BLtree_scaledF_entire <- predvars2deciles(df = preds_trees %>% select(AngioTreeCover_prop, percBLTree_scaled_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na()
                                   , response_vars = c("AngioTreeCover_prop", "percBLTree_scaled_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_scaledF_entire <- decile_dotplot_pq(df = deciles_BLtree_scaledF_entire, response= c("AngioTreeCover_prop", "percBLTree_scaled_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Proportion of trees that are Broad-leaved (scaled to sum to 1 w/ BL tree) - forest model in CONUS")

quantPlot_BLtree_scaledF_entire <- add_dotplot_inset(quantPlot_BLtree_scaledF_entire, df = deciles_BLtree_scaledF_entire, response = c("AngioTreeCover_prop", "percBLTree_scaled_F"), dfRaw = preds_trees %>% select(AngioTreeCover_prop, percBLTree_scaled_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

 
## prediction of absolute cover of BL tree in forest
# forest model predicting in forest (EPA-defined)
deciles_BLtree_absoluteF_F <- predvars2deciles(df = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(BLTreeCover_abs, absBLTree_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na() 
                                   , response_vars = c("BLTreeCover_abs", "absBLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_absoluteF_F <- decile_dotplot_pq(df = deciles_BLtree_absoluteF_F, response= c("BLTreeCover_abs", "absBLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute cover of Broad-leaved trees - forest model in forest")

quantPlot_BLtree_absoluteF_F <- add_dotplot_inset(quantPlot_BLtree_absoluteF_F, df = deciles_BLtree_absoluteF_F, response = c("BLTreeCover_abs", "absBLTree_F"), dfRaw = preds_trees %>% filter(newRegion != "dryShrubGrass") %>% select(BLTreeCover_abs, absBLTree_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

# forest model predicting in CONUS
deciles_BLtree_absoluteF_entire <- predvars2deciles(df = preds_trees %>% select(BLTreeCover_abs, absBLTree_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na()
                                   , response_vars = c("BLTreeCover_abs", "absBLTree_F"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse"), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_absoluteF_entire <- decile_dotplot_pq(df = deciles_BLtree_absoluteF_entire, response= c("BLTreeCover_abs", "absBLTree_F"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute cover of Broad-leaved trees - forest model in CONUS")

quantPlot_BLtree_absoluteF_entire <- add_dotplot_inset(quantPlot_BLtree_absoluteF_entire, df = deciles_BLtree_absoluteF_entire, response = c("BLTreeCover_abs", "absBLTree_F"), dfRaw = preds_trees %>% select(BLTreeCover_abs, absBLTree_F, "tmean", "prcp_dry", "prcpTempCorr", "sand", "carbon", "AWHC", "prcpTempCorr_anom", "isothermality", "isothermality_anom", "coarse") %>% drop_na(), add_smooth = TRUE, deciles = FALSE)


## prediction of absolute cover of BL tree, scaled by ecoregion probability to be CONUS-wide 
deciles_BLtree_absolute_CONUS <- predvars2deciles(df = preds_trees %>% select(BLTreeCover_abs, absBLTree_CONUS, "tmean",  "prcp" , "prcpTempCorr" , "sand"  ,   "coarse"   , "AWHC",   "prcp_seasonality_anom" ,"isothermality_anom" ,  "isothermality"  ,       "prcp_seasonality",  "prcpTempCorr_anom" ,    "prcp_dry" ,   "carbon" 
                                                                              
                                                                              ) %>% drop_na()
                                   , response_vars = c("BLTreeCover_abs", "absBLTree_CONUS"), # name of observations, followed by name of predictions
                                   pred_vars = c("tmean",  "prcp" , "prcpTempCorr" , "sand"  ,   "coarse"   , "AWHC",   "prcp_seasonality_anom" ,"isothermality_anom" ,  "isothermality"  ,       "prcp_seasonality",  "prcpTempCorr_anom" ,    "prcp_dry" ,   "carbon" ), # for predictors, combine list of predictors for all models (total tree for F and F, and perc BL tree for F and F)
cut_points = seq(0, 1, 0.01)
                                     )
quantPlot_BLtree_absolute_CONUS <- decile_dotplot_pq(df = deciles_BLtree_absolute_CONUS, response= c("BLTreeCover_abs", "absBLTree_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute cover of Broad-leaved trees - forest model in CONUS")

quantPlot_BLtree_absolute_CONUS <- add_dotplot_inset(quantPlot_BLtree_absolute_CONUS, df = deciles_BLtree_absolute_CONUS, response = c("BLTreeCover_abs", "absBLTree_CONUS"), dfRaw = preds_trees %>% select(BLTreeCover_abs, absBLTree_CONUS, "tmean",  "prcp" , "prcpTempCorr" , "sand"  ,   "coarse"   , "AWHC",   "prcp_seasonality_anom" ,"isothermality_anom" ,  "isothermality"  ,       "prcp_seasonality",  "prcpTempCorr_anom" ,    "prcp_dry" ,   "carbon" ) %>% drop_na(), add_smooth = TRUE, deciles = FALSE)

#### Get quantile plots in order for google doc ####
### Total herb
## Forest model
#Predictions of absolute total tree cover (absTotalTree_F) in forest (EPA ecoregion-defined) using the forest model
quantPlot_totalHerb_F_F
#Predictions of absolute total tree cover (absTotalTree_F) in CONUS using the forest model
quantPlot_totalHerb_F_entire
## Grass/shrub model
# Predictions of absolute total tree cover (absTotalTree_GS) in grass/shrub (EPA ecoregion-defined) using the grass/shrub model
quantPlot_totalHerb_GS_GS
# Predictions of absolute total tree cover (absTotalTree_GS) in CONUS using the grass/shrub model
quantPlot_totalHerb_GS_entire

### total trees
## Forest model
#Predictions of absolute total tree cover (absTotalTree_F) in forest (EPA ecoregion-defined) using the forest model 
quantPlot_totalTree_F_F
#Predictions of absolute total tree cover (absTotalTree_F) in forest (model-defined) using the forest model 
quantPlot_totalTree_F_Fprob
#Predictions of absolute total tree cover (absTotalTree_F) in CONUS using the forest model
quantPlot_totalTree_F_entire
## Grass/shrub model
# Predictions of absolute total tree cover (absTotalTree_GS) in grass/shrub (EPA ecoregion-defined) using the grass/shrub model
# quantPlot_totalTree_GS_GS
# # Predictions of absolute total tree cover (absTotalTree_GS) in grass/shrub (model-defined) using the grass/shrub model
# quantPlot_totalTree_GS_GSprob
# # Predictions of absolute total tree cover (absTotalTree_GS) in CONUS using the grass/shrub model
# quantPlot_totalTree_GS_entire

### NL trees
##Forest model
#Predictions of proportion of trees that are needle-leaved (percNLTree_F) in forest (EPA ecoregion-defined) using the forest model
quantPlot_NLtree_RawF_F
#Predictions of proportion of trees that are needle-leaved (percNLTree_F) in forest (model-defined) using the forest model 
quantPlot_NLtree_RawF_Fprob
#Predictions of proportion of trees that are needle-leaved (percNLTree_F) in CONUS using the forest model
quantPlot_NLtree_RawF_entire

##Grass/shrub model
#Predictions of proportion of trees that are needle-leaved (percNLTree_GS) in grass/shrub (EPA ecoregion-defined) using the grass/shrub model
# quantPlot_NLtree_RawGS_GS
# #Predictions of proportion of trees that are needle-leaved (percNLTree_GS) in grass/shrub (model-defined) using the grass/shrub model
# quantPlot_NLtree_RawGS_GSprob
# #Predictions of proportion of trees that are needle-leaved (percNLTree_GS) in CONUS using the grass/shrub model
# quantPlot_NLtree_RawGS_entire

## scaling
#Predictions of proportion of trees that are needle-leaved in Grass/Shrub, scaled w/ broad-leaved tree in GS to sum to 1 (percNLTree_scaled_GS = percNLTree_GS/(percBLTree_GS + percNLTree_GS))
# quantPlot_NLtree_scaledGS_GS
# quantPlot_NLtree_scaledGS_entire
# Predictions of proportion of trees that are needle-leaved in Forest, scaled w/ broad-leaved tree in F to sum to 1 (percNLTree_scaled_F = percNLTree_F/(percBLTree_F + percNLTree))
quantPlot_NLtree_scaledF_F
quantPlot_NLtree_scaledF_entire
# Predictions of absolute cover of needle-leaved in Grass/Shrub (absNLTree_GS = (percNLTree_scaled_GS * absTotalTree_GS) )
# quantPlot_NLtree_absoluteGS_GS
# quantPlot_NLtree_absoluteGS_entire
# Predictions of  absolute cover of trees that are needle-leaved in Forest (absNLTree_F =(percNLTree_scaled_F * absTotalTree_F))
quantPlot_NLtree_absoluteF_F
quantPlot_NLtree_absoluteF_entire
# Scale the predictions of proportion of trees that are needle-leaved CONUS-wide (Equation: absNLTree_CONUS = (prob_Forest * absNLTree_F)  + (prob_grassShrub * absNLTree_GS))
quantPlot_NLtree_absolute_CONUS

### BL trees
##Forest model
#Predictions of proportion of trees that are broad-leaved (percBLTree_F) in forest (EPA ecoregion-defined) using the forest model
quantPlot_BLtree_RawF_F
#Predictions of proportion of trees that are broad-leaved (percBLTree_F) in forest (model-defined) using the forest model 
quantPlot_BLtree_RawF_Fprob
#Predictions of proportion of trees that are broad-leaved (percBLTree_F) in CONUS using the forest model
quantPlot_BLtree_RawF_entire

##Grass/shrub model
#Predictions of proportion of trees that are broad-leaved (percBLTree_GS) in grass/shrub (EPA ecoregion-defined) using the grass/shrub model
# quantPlot_BLtree_RawGS_GS
# #Predictions of proportion of trees that are broad-leaved (percBLTree_GS) in grass/shrub (model-defined) using the grass/shrub model
# quantPlot_BLtree_RawGS_GSprob
# #Predictions of proportion of trees that are broad-leaved (percBLTree_GS) in CONUS using the grass/shrub model
# quantPlot_BLtree_RawGS_entire

## scaling
#Predictions of proportion of trees that are broad-leaved in Grass/Shrub, scaled w/ broad-leaved tree in GS to sum to 1 (percBLTree_scaled_GS = percBLTree_GS/(percBLTree_GS + percBLTree_GS))
# quantPlot_BLtree_scaledGS_GS
# quantPlot_BLtree_scaledGS_entire
# Predictions of proportion of trees that are broad-leaved in Forest, scaled w/ broad-leaved tree in F to sum to 1 (percBLTree_scaled_F = percBLTree_F/(percBLTree_F + percBLTree))
quantPlot_BLtree_scaledF_F
quantPlot_BLtree_scaledF_entire
# Predictions of absolute cover of broad-leaved in Grass/Shrub (absBLTree_GS = (percBLTree_scaled_GS * absTotalTree_GS) )
# quantPlot_BLtree_absoluteGS_GS
# quantPlot_BLtree_absoluteGS_entire
# Predictions of  absolute cover of trees that are broad-leaved in Forest (absBLTree_F =(percBLTree_scaled_F * absTotalTree_F))
quantPlot_BLtree_absoluteF_F
quantPlot_BLtree_absoluteF_entire
# Scale the predictions of proportion of trees that are broad-leaved CONUS-wide (Equation: absBLTree_CONUS = (prob_Forest * absBLTree_F)  + (prob_grassShrub * absBLTree_GS))
quantPlot_BLtree_absolute_CONUS
#### Compare raw data for absolute needle leaved tree cover ####
# observations in grass/shrub
ggplot(preds_trees %>% filter(newRegion == "dryShrubGrass")) + 
  geom_density(aes(TotalTreeCover)) + # observed total tree cover
  geom_density(aes(ConifTreeCover_prop), col = "orange") + # % of needle-leaved tree 
  geom_density(aes(NLTreeCover_abs), col = "red") + # observed absolute NL tree cover
  #geom_density(aes(TotalTreeCover * ConifTreeCover_prop), col = "grey") + 
  labs(title = "Observed needle-leaved tree cover = Grass/Shrub", subtitle = "black = total tree cover \norange = % of tree cover that is needle-leaved \nred = absolute needle-leaved tree cover")
# observations in forest
ggplot(preds_trees %>% filter(newRegion != "dryShrubGrass")) + 
  geom_density(aes(TotalTreeCover)) + # observed total tree cover
  geom_density(aes(ConifTreeCover_prop), col = "orange") + # % of needle-leaved tree 
  geom_density(aes(NLTreeCover_abs), col = "red") + # observed absolute NL tree cover
  #geom_density(aes(TotalTreeCover * ConifTreeCover_prop), col = "grey") + 
  labs(title = "Observed needle-leaved tree cover = Forest", subtitle = "black = total tree cover \norange = % of tree cover that is needle-leaved \nred = absolute needle-leaved tree cover")
  
#grass/shrub
ggplot(preds_trees) + 
    geom_density(aes(absTotalTree_GS), col = "black") +# modeled total tree cover
  geom_density(aes(percNLTree_scaled_GS), col = "orange") + # % of needle-leaved tree (scaled to CONUS-wide)
    geom_density(aes(absNLTree_GS), col = "red") + # modeled absolute NL tree cover
  #geom_density(aes(absTotalTree_CONUS * percNLTree_scaled_CONUS), col = "grey") + 
  labs(title = "Modeled needle-leaved tree cover - Grass/Shrub Ecoregion", subtitle = "black = total tree cover \norange = % of tree cover that is needle-leaved  \nred = absolute needle-leaved tree cover ")
#forest
ggplot(preds_trees) + 
    geom_density(aes(absTotalTree_F), col = "black") +# modeled total tree cover
  geom_density(aes(percNLTree_scaled_F), col = "orange") + # % of needle-leaved tree (scaled to CONUS-wide)
    geom_density(aes(absNLTree_F), col = "red") + # modeled absolute NL tree cover
  #geom_density(aes(absTotalTree_CONUS * percNLTree_scaled_CONUS), col = "grey") + 
  labs(title = "Modeled needle-leaved tree cover - Forest Ecoregion", subtitle = "black = total tree cover \norange = % of tree cover that is needle-leaved  \nred = absolute needle-leaved tree cover ")


## comparing distributions of raw values for NL trees by hand 
# scaled percentage of NL tree in forest - grass/shrub
violin_GS_obs <- preds_trees %>% 
  filter(prob_Forest < 0.5) %>% 
  select(ConifTreeCover_prop, percNLTree_scaled_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% 
  drop_na() %>% 
  pivot_longer(cols = c(tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse)) %>%
  ggplot() +
  facet_wrap(~name, scales = "free") +
  geom_point(aes(x = value, y = ConifTreeCover_prop), alpha = .1, col = "grey") + # observed
 # geom_point(aes(x = value, y = percNLTree_scaled_GS), alpha = .1, col = "lightblue") + # predicted
  geom_violin(aes(x = value, y = ConifTreeCover_prop, group = cut_number(value,n = 10)), col = "black") #+ # observed 
  #geom_violin(aes(x = jitter(value), y = percNLTree_scaled_GS, group = cut_number(value,n = 10)), col = "blue") # predicted
violin_GS_preds <- preds_trees %>% 
  filter(prob_Forest < 0.5) %>% 
  select(ConifTreeCover_prop, percNLTree_scaled_GS, tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse) %>% 
  drop_na() %>% 
  pivot_longer(cols = c(tmean, prcp, AWHC, prcpTempCorr, isothermality, prcpTempCorr_anom, coarse)) %>%
  ggplot() +
  facet_wrap(~name, scales = "free") +
  #geom_point(aes(x = value, y = ConifTreeCover_prop), alpha = .1, col = "grey") + # observed
  geom_point(aes(x = value, y = percNLTree_scaled_GS), alpha = .1, col = "lightblue") + # predicted
  #geom_violin(aes(x = value, y = ConifTreeCover_prop, group = cut_number(value,n = 10)), col = "black") #+ # observed 
  geom_violin(aes(x = jitter(value), y = percNLTree_scaled_GS, group = cut_number(value,n = 10)), col = "blue") # predicted
ggarrange(violin_GS_obs, violin_GS_preds, ncol = 1)


### locations where we have cover for total trees vs locations where we have values for cover of NL or BL trees
totTreeLocations <- preds_quantile %>% 
  select(TotalTreeCover, x, y) %>% 
  drop_na() %>% 
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "TotalTreeCover", 
                   fun = mean, na.rm = TRUE) 
ggarrange(
  ggplot() +
  geom_spatraster(data = totTreeLocations) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Observations of total cover (all available points)")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])
  ,
  ggplot() + 
    geom_density(data = data.frame("TotalTreeCover" = as.numeric(na.omit(terra::values(totTreeLocations$mean)))), aes(x = TotalTreeCover)),
ncol = 1
)

totTreeLocations_trim <- preds_quantile %>% 
  select(TotalTreeCover, AngioTreeCover_prop, ConifTreeCover_prop, x, y) %>% 
  drop_na() %>% 
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "TotalTreeCover", 
                   fun = mean, na.rm = TRUE) 
ggarrange(
  ggplot() +
  geom_spatraster(data = totTreeLocations_trim) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Observations of total cover (locations w/ NL and BL cover)")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])
  ,
  ggplot() + 
    geom_density(data = data.frame("TotalTreeCover" = as.numeric(na.omit(terra::values(totTreeLocations_trim$mean)))), aes(x = TotalTreeCover)),
ncol = 1
)
# 
### locations where we have cover for total herbaceous vs locations where we have values for cover of C3/C4/forbs
totHerbLocations <- preds_quantile %>% 
  select(TotalHerbaceousCover, x, y) %>% 
  drop_na() %>% 
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "TotalHerbaceousCover", 
                   fun = mean, na.rm = TRUE) 
ggarrange(
  ggplot() +
  geom_spatraster(data = totHerbLocations) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Observations of total herbaceous cover (all available points)")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])
  ,
  ggplot() + 
    geom_density(data = data.frame("TotalHerbaceousCover" = as.numeric(na.omit(terra::values(totHerbLocations$mean)))), aes(x = TotalHerbaceousCover)),
ncol = 1
)

totHerbLocations_trim <- preds_quantile %>% 
  select(TotalHerbaceousCover, C3GramCover_prop, C4GramCover_prop, ForbCover_prop, x, y) %>% 
  drop_na() %>% 
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "TotalHerbaceousCover", 
                   fun = mean, na.rm = TRUE) 
ggarrange(
  ggplot() +
  geom_spatraster(data = totHerbLocations_trim) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Observations of total herbaceous cover (locations w/ C3/C4/forb cover)")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])
  ,
  ggplot() + 
    geom_density(data = data.frame("TotalHerbaceousCover" = as.numeric(na.omit(terra::values(totHerbLocations_trim$mean)))), aes(x = TotalHerbaceousCover)),
ncol = 1
)

# total trees trees 
# absolute cover
quantPlot_totalTrees <- decile_dotplot_pq(df = deciles_totTree, response= c("TotalTreeCover", "absTotalTree_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Tree Cover: Quantile Plot")

quantPlot_totalTrees <- add_dotplot_inset(quantPlot_totalTrees, df = deciles_totTree, response = c("TotalTreeCover", "absTotalTree_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# needle-leaved trees 
# absolute cover
quantPlot_NLtreeAbs <- decile_dotplot_pq(df = deciles_NLtree, response= c("NLTreeCover_abs", "absNLTree_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Needle-Leaved Tree Cover: Quantile Plot")

quantPlot_NLtreeAbs <- add_dotplot_inset(quantPlot_NLtreeAbs, df = deciles_NLtree, response = c("NLTreeCover_abs", "absNLTree_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# broad-leaved trees 
# absolute cover
quantPlot_BLtreeAbs <- decile_dotplot_pq(df = deciles_BLtree, response= c("BLTreeCover_abs", "absBLTree_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Broad-Leaved Tree Cover: Quantile Plot")

quantPlot_BLtreeAbs <- add_dotplot_inset(quantPlot_BLtreeAbs, df = deciles_BLtree, response = c("BLTreeCover_abs", "absBLTree_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# shrubs 
quantPlot_shrubs <- decile_dotplot_pq(df = deciles_shrub, response= c("ShrubCover", "absShrub_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Shrub Cover: Quantile Plot")

quantPlot_shrubs <- add_dotplot_inset(quantPlot_shrubs, df = deciles_shrub, response = c("ShrubCover", "absShrub_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# bare ground
quantPlot_bareGround <- decile_dotplot_pq(df = deciles_bareGround, response= c("BareGroundCover", "absBareGround_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Bare Ground Cover: Quantile Plot")

quantPlot_bareGround <- add_dotplot_inset(quantPlot_bareGround, df = deciles_bareGround, response = c("BareGroundCover", "absBareGround_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# total herbaceous
quantPlot_totHerb <- decile_dotplot_pq(df = deciles_totalHerb, response= c("TotalHerbaceousCover", "absTotalHerb_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Herbaceous Cover: Quantile Plot")

quantPlot_totHerb <- add_dotplot_inset(quantPlot_totHerb, df = deciles_totalHerb, response = c("TotalHerbaceousCover", "absTotalHerb_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# forbs
quantPlot_forbs <- decile_dotplot_pq(df = deciles_forb, response= c("ForbCover_abs", "absForb_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Forb Cover: Quantile Plot")

quantPlot_forbs <- add_dotplot_inset(quantPlot_forbs, df = deciles_forb, response = c("ForbCover_abs", "absForb_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# C3 grass
quantPlot_c3grass <- decile_dotplot_pq(df = deciles_C3grass, response= c("C3GramCover_abs", "absC3grass_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute C3 Grass Cover: Quantile Plot")

quantPlot_c3grass <- add_dotplot_inset(quantPlot_c3grass, df = deciles_C3grass, response = c("C3GramCover_abs", "absC3grass_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

#C4 grass 
quantPlot_c4grass <- decile_dotplot_pq(df = deciles_C4grass, response= c("C4GramCover_abs", "absC4grass_CONUS"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute C4 Grass Cover: Quantile Plot")

quantPlot_c4grass <- add_dotplot_inset(quantPlot_c4grass, df = deciles_C4grass, response = c("C4GramCover_abs", "absC4grass_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

## now, display the figures
quantPlot_totalTrees

quantPlot_totHerb

quantPlot_shrubs

quantPlot_bareGround

quantPlot_NLtreeAbs

quantPlot_BLtreeAbs

quantPlot_c3grass

quantPlot_c4grass

quantPlot_forbs

Below, we show the RMSE and bias of the final, scaled estimates of absolute cover for each functional type in comparison to the observed data.

long_absPreds <- preds_quantile %>% 
  mutate(uniqueID = 1:nrow(preds_quantile)) %>% 
  select(uniqueID, absTotalHerb_CONUS, absTotalTree_CONUS, absShrub_CONUS, absBareGround_CONUS, 
                        absNLTree_CONUS, absBLTree_CONUS, absC3grass_CONUS, absC4grass_CONUS, absForb_CONUS) %>% 
  rename(TotalHerb = absTotalHerb_CONUS, 
         TotalTree = absTotalTree_CONUS, 
         Shrub = absShrub_CONUS, 
         BareGround = absBareGround_CONUS, 
         NLTree = absNLTree_CONUS, 
         BLTree = absBLTree_CONUS, 
         C3gram = absC3grass_CONUS, 
         C4gram = absC4grass_CONUS, 
         Forb = absForb_CONUS) %>% 
  pivot_longer(cols = c(TotalHerb, TotalTree, Shrub, BareGround, NLTree, BLTree, C3gram, C4gram, Forb), 
                names_to = "names", 
               values_to = "preds_values"
               ) %>% 
  select(uniqueID, names, preds_values)

long_absObs <- preds_quantile %>% 
  mutate(uniqueID = 1:nrow(preds_quantile)) %>% 
  select(uniqueID, TotalHerbaceousCover, TotalTreeCover, ShrubCover, BareGroundCover,
                        NLTreeCover_abs, BLTreeCover_abs, C3GramCover_abs, C4GramCover_abs, ForbCover_abs) %>% 
  rename(TotalHerb = TotalHerbaceousCover, 
         TotalTree = TotalTreeCover, 
         Shrub = ShrubCover, 
         BareGround = BareGroundCover, 
         NLTree = NLTreeCover_abs, 
         BLTree = BLTreeCover_abs, 
         C3gram = C3GramCover_abs, 
         C4gram = C4GramCover_abs, 
         Forb = ForbCover_abs) %>% 
  pivot_longer(cols = c(TotalHerb, TotalTree, Shrub, BareGround, NLTree, BLTree, C3gram, C4gram, Forb), 
                names_to = "names", 
               values_to = "obs_values"
               ) %>% 
  select(uniqueID, names, obs_values)

longObsPreds <- long_absObs %>% 
  left_join(long_absPreds)

## calculate bias
(summaryTableAbs <- longObsPreds %>% 
  mutate(resids = obs_values - preds_values) %>% 
  group_by(names) %>% 
  dplyr::summarise(bias = mean(resids, na.rm = TRUE),
            RMSE = round(sqrt(mean(resids^2, na.rm = TRUE)),4)
            ) %>% 
  slice(8, 9, 7, 2, 6, 1, 3, 4, 5) %>% kable(
  format = "html", 
  col.names = c("Model Name", "Model bias mean(obs. - pred.)", "model RMSE")
))

Model Name	Model bias mean(obs. - pred.)	model RMSE
TotalHerb	0.0006038	0.1443
TotalTree	0.0214052	0.0974
Shrub	0.0000854	0.0867
BareGround	-0.0000675	0.0918
NLTree	0.1066927	0.2295
BLTree	0.0336315	0.1674
C3gram	0.0240471	0.1557
C4gram	-0.0024318	0.0745
Forb	-0.0176334	0.1250

Finally, we show how these model predictions made using contemporary climate and weather data change over time from 2010 to 2020.

# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDat_time_unscaled <- climDat_time
preds_byHand_time <- climDat_time_unscaled %>% 
  mutate(pred =  as.numeric(1/ (1 + exp(-(-0.67478165 +  
                                            1.71345065*((prcp - 483.34337211) / 321.52170566) + 
                                            -0.01633932*((prcp_seasonality - 0.97338094) / 0.23365057) + 
                                            -0.91285243*((prcpTempCorr - 0.01681097) / 0.39882063) +  
                                            0.53377439*((annWetDegDays - 1808.61043539) / 1083.25298131) + 
                                            -0.60523936*((annWatDef_95 - 164.68611646) / 104.14892564) +  
                                            0.20578879*((coarse - 9.66120017) / 10.16883531) + 
                                            -0.61568228*((AWHC - 14.61190409) / 5.65765542) + 
                                            -0.31818631*I(((tmean - 11.0291275) / 4.89048541)^2) + 
                                            0.20104126*I(((prcp_seasonality - 0.97338094) / 0.23365057)^2) + 
                                            -0.89055646*I(((prcpTempCorr - 0.01681097) / 0.39882063)^2) +  
                                            0.14226917*I(((isothermality - 37.76731082) / 5.2256935)^2) + 
                                            -0.11810211*I(((clay - 19.59171509) / 8.51944876)^2) +  
                                            0.11278028*I(((AWHC - 14.61190409) / 5.65765542)^2) +  
                                            0.91656456*((annWetDegDays - 1808.61043539) / 1083.25298131)*((annWatDef_95 - 164.68611646) / 104.14892564) + 
                                            -0.29650844*((prcp_dry - 3.18622575) / 5.36400552)*((isothermality - 37.76731082) / 5.2256935) + 
                                            -0.02061048*((prcpTempCorr - 0.01681097) / 0.39882063)*((isothermality - 37.76731082) / 5.2256935) +
                                            -0.28504255*((isothermality - 37.76731082) / 5.2256935)*((tmean - 11.0291275) / 4.89048541) + 
                                            0.50549093*((prcpTempCorr - 0.01681097) / 0.39882063)*((prcp_dry - 3.18622575) / 5.36400552) +
                                            -0.25285633*((prcpTempCorr - 0.01681097) / 0.39882063)*((tmean - 11.0291275) / 4.89048541) + 
                                            -0.16290284*((AWHC - 14.61190409) / 5.65765542)*((clay - 19.59171509) / 8.51944876) +  
                                            0.09163249*((coarse - 9.66120017) / 10.16883531)*((clay - 19.59171509) / 8.51944876))))))
# convert the continuous probability predictions to binary using previously-identified threshold 
preds_byHand_time <- preds_byHand_time %>% 
  mutate(pred_binary = pred,
    pred_binary = replace(pred_binary, pred_binary > 0.3819095, 1),
         pred_binary = replace(pred_binary, pred_binary <= 0.3819095, 0))
# interpolate continuous predictions
preds_byHand_time$preds_interp <- NA
preds_byHand_time[!is.na(preds_byHand_time$pred), "preds_interp"] <- approxfun_pred(preds_byHand_time[ !is.na(preds_byHand_time$pred), "pred"])

# names(climDat_time_unscaled)[c(5, 6, 7, 13, 10, 12, 99, 101, 102, 103)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
# 
# # predict with contemporary climate data 
# preds_byHand_time <- climDat_time_unscaled %>% 
#   mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM +  0.2456*T_coldestMonth_meanAnnAvg_CLIM +  0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM +  0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth +  0.0335*avgSandPerc_acrossDepth +  0.0310*avgCoarsePerc_acrossDepth +  0.2726*avgOrganicCarbonPerc_0_3cm))))
# #predictionsModel <- predict(object = forestMod, type = "response", newdata = climDatPred_unscaled)

## with contemporary climate data
totalHerb_F_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsOverTime <-
    makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalHerb_GS) %>%
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsOverTime<- climDat_time_Pred %>% 
  mutate(absTotalTree_GS = 0)
  # makePredictions(predictionDF = climDat_time_Pred,
  #                                           modelObject = totalTree_GS) %>% 
  # rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsOverTime<- climDat_time_Pred %>% 
  mutate(percBLTree_GS = 0)
  # makePredictions(predictionDF = climDat_time_Pred,
  #                                           modelObject = BLtree_GS) %>% 
  # rename(percBLTree_GS = "modelPreds")

NLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsOverTime<- climDat_time_Pred %>% 
  mutate(percNLTree_GS = 0)
  # makePredictions(predictionDF = climDat_time_Pred,
  #                                           modelObject = NLtree_GS) %>% 
  # rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

predsOverTime <- totalHerb_F_predsOverTime  %>% 
  cbind(totalHerb_GS_predsOverTime %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsOverTime %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsOverTime %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsOverTime %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsOverTime %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsOverTime %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsOverTime %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsOverTime %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsOverTime %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsOverTime %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsOverTime %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsOverTime %>% select(percForb_CONUS)) %>% 
  cbind(preds_byHand_time %>% select(preds_interp)) %>% 
  rename(prob_Forest = "preds_interp") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

## now, scale predictions according to ecoregion percentage
predsOverTime <- predsOverTime %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
         percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
         percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
         )
## then, scale the percentages so they sum to 1
predsOverTime <- predsOverTime %>% 
  mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

## convert percentages into absolute cover for level 2 cover classes
predsOverTime <- predsOverTime %>% 
    mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS), 
         absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
  )

# change spatialIDs to discrete numbers, since the long form is causing some sort of floating point issue? 
spatID_lu <- data.frame("spatialID" = unique(predsOverTime$spatialID), 
                        "uniqueID" = 1:length(predsOverTime$spatialID))
predsOverTime <- predsOverTime %>% 
  left_join(spatID_lu)

# make into a long format
predsOverTime_long <- predsOverTime %>% 
  select(year,x, y, spatialID, uniqueID, absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS) %>% 
  pivot_longer(cols = c(absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS)) %>% 
  mutate(year = as.integer(year))

# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>% 
  group_by(uniqueID, name) %>% 
  summarize(cv = sd(value)/mean(value))
## add x,y 
predsOverTime_cv <- predsOverTime_cv %>% 
  left_join(predsOverTime %>% select(uniqueID, x, y))

# visualize change over time for a selection of points
coverOverTimeForSelectLocations <- predsOverTime_long %>% 
  filter(uniqueID %in% 1:50) %>% 
  filter(name %in% c("absShrub_CONUS", "absBareGround_CONUS",  "absNLTree_CONUS", 
                     "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")) %>% 
  ggplot() + 
  facet_wrap(~uniqueID) + 
  geom_area(aes(x = year, y = value, fill = name)) + 
  geom_line(data = predsOverTime_long %>% 
              filter(uniqueID %in% 1:50) %>% 
              filter(name %in% c("absTotalHerb_CONUS", "absTotalTree_CONUS")),
            aes(x = year, y = value, linetype = name), col = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90)) + 
  labs(title = "Change in modeled absolute cover by functional type from \n 2010 to 2020 for a subset of locations")
mapOfSelectLocations <- predsOverTime_long %>% 
    select(x, y, uniqueID) %>% 
    filter(uniqueID %in% 1:50) %>% 
    unique() %>% 
    ggplot() + 
    geom_sf(data = cropped_states, aes()) + 
    geom_point(aes(x, y)) + 
    geom_label_repel(aes(x, y, label = uniqueID)) + 
    labs(title = "Locations in above plot")
    

ggarrange(coverOverTimeForSelectLocations, mapOfSelectLocations, ncol = 1, heights = c(1,.5))

predsOverTime_cv %>%  
  ggplot() + 
  geom_boxplot(aes(x = name, y = cv, col = name)) +
  theme_minimal() + 
  labs(title = "Coefficients of variation of model-predicted absolute \n cover at 1000 random locations for each functional group from 2010-2020",
        x = "functional group",
       y = "coefficient of variation") +
  theme(axis.text.x = element_text(angle = 90))

# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>% 
 # rename(geometry = Shape)
 library(hexbin)
ggplot() + 
  facet_wrap(~name) + 
  geom_sf(data = cropped_states, aes()) + 
  stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
  geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) + 
  guides(fill = guide_legend(title = "CV")) + 
  theme_minimal() + 
  labs(title = "The spatial distribution of the coefficients of variation of \nmodel-predicted cover from 2010-2020 for 1000 random locations",
       subtitle = "Note that values are aggregated accross space to enhance readability")

Now, we generate CONUS-wide model predictions of relative cover

We use the model predictions of absolute cover we generated above, which are CONUS-wide and have been scaled according to predicted ecoregions, to calculate relative cover for each functional group. We consider the overstory (trees) and understory (all other functional types) separately in our relativization scheme, such that the values of relative cover are identical to the values of absolute cover for total trees, needle-leaved trees, and broad-leaved trees. For the understory (comprised of all herbaceous functional types, shrubs, and bare ground), total relative cover must sum to 1, so absolute cover for these functional groups is scaled accordingly.

## for contemporary data
contempPreds <- contempPreds %>% 
  mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                           absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
  mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                           absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
  mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                           absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

Below are maps of predicted relative cover for each functional group we modeled. These maps show model-predicted relative cover using contemporary climate and weather data, residuals when comparing these contemporary predictions to observations, and changes in model predictions of cover when predictions are made using simulated future climate and weather data from two GCMs. This figure also includes a map of the total understory absolute cover, which was the scaling factor used to calculate relative cover for each understory group. Note that relative cover maps for total trees, needle-leaved trees, and broad-leaved trees are not shown, since, for those functional groups, relative cover is identical to absolute cover.

# list of absolute cover names
responseNames <- c(#"relCoverB_totalTree", 
  "relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
                   #"relCoverB_NLTree", "relCoverB_BLTree", 
  "relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
                   "total_NonTree_AbsCover_CONUS")

# prepare the observed data for raster calculations (calculate relative observed cover)
obsDat <- obsDat %>% 
  mutate(
       total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + absShrub_CONUS + absBareGround_CONUS) # for observations
       ) %>% 
  mutate(
  # making observations relative
         relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS,
         )
relCover_B_Maps_contemp <- lapply(responseNames, 
                                      FUN = function(x) {
                                      
                                       ## contemporary absolute cover
                                        # rasterize response 
                                        resp_rast <- contempPreds %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_2 <- resp_rast %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        
                                        ## now, rasterize the observed data and calculate residualks 
                                         plotObservations <- obsDat %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs("EPSG:4326")) %>% 
                                          terra::project(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE) 
                                         plotObservations <- plotObservations %>%  
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate residuals (observed - predicted)
                                        resids_rast <- plotObservations - resp_rast_2
                                        ## aggregate the residuals raster to make it easier to see 
                                        resids_rastCoarse <- resids_rast %>% 
                                          aggregate(fact = 2, fun = "mean", na.rm = TRUE)
                                        # map the residuals
                                        
                                        mapResids <- ggplot() +
                                          geom_spatraster(data = resids_rastCoarse) + 
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% 
                                                    st_crop(st_bbox(resids_rast)),fill=NA ) +
                                          labs(title = paste0("Residuals (observations - predicted absolute cover) for ",x)) +
                                          scale_fill_gradient2(low = "red",
                                                               mid = "white" ,
                                                               high = "blue" , 
                                                               midpoint = 0,   na.value = "lightgrey",
                                                               limits = c(-1,1)
                                          )  + 
                                          xlim(st_bbox(resids_rast)[c(1,3)]) + 
                                          ylim(st_bbox(resids_rast)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram of residuals
                                        histResids <- ggplot() + 
                                          geom_density(aes(terra::values(resids_rast$mean, na.rm = TRUE, mat = FALSE)), 
                                                       col = "black", fill = "darkgrey") + 
                                          geom_vline(aes(xintercept = 0), linetype = 2) +
                                          xlab(paste0("Residuals (obs - pred) ",x)) +
                                          theme_classic()
                                        
                                         ## future model 1 absolute cover
                                        # rasterize response 
                                        resp_rast_future1 <- predsFuture1 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future1_2 <- resp_rast_future1 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
                                        
                                        ## future model 2 absolute cover
                                        # rasterize response 
                                        resp_rast_future2 <- predsFuture2 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future2_2 <- resp_rast_future2 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
                                        
                                        
                                        if (x == "total_NonTree_AbsCover_CONUS") {
                                          # make a map of predictions w/ current data 
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
                                          scale_fill_gradient2(low = "lightblue",
                                                               #mid = "wheat" ,
                                                               high = "darkblue" , 
                                                               #midpoint = 0, 
                                                               limits = c(0,max(terra::values(resp_rast_2$mean), na.rm= TRUE)),  
                                                              na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
                                          theme_classic()
                                        
                                         # make a map of future preds w/ model 1
                                         
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of Total Absolute Cover (non-tree) \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-2,2),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        # make a map of future preds w/ model 2
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of Cover (non-tree), \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-2,2),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        
                                        } else {
                                          # make a map of predictions w/ current data 
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of relative cover of ",x,
                                                                        ", \n when all functional groups except trees sum to 1")) +
                                          scale_fill_gradient2(low = "brown",
                                                               mid = "wheat" ,
                                                               high = "darkgreen" , 
                                                               midpoint = 0, limits = c(0,1),  
                                                              na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x)) +
                                          theme_classic()
                                        
                                         # make a map of future preds w/ model 1
                                         
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of relative cover of ",x,", \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        # make a map of future preds w/ model 2
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of relative cover of ",
                                                                        x,
                                                                        ", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        }
                                        
                                        
                                        return(list("rast" = resp_rast_2, 
                                                    "map" = resp_map, 
                                                    "hist" = resp_hist,
                                                    "rast_resids" = resids_rast,
                                                    "map_resids" = mapResids,
                                                    "hist_resids" = histResids,
                                                    "rast_future1" = resp_rast_future1_2,
                                                    "map_future1" = resp_map_future1, 
                                                    "hist_future1" = resp_hist_future1,
                                                    "rast_future2" = resp_rast_future2_2,
                                                    "map_future2" = resp_map_future2, 
                                                    "hist_future2" = resp_hist_future2))
                                        })


names(relCover_B_Maps_contemp) <- c(#"relCoverB_totalTree", 
  "relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
                   #"relCoverB_NLTree", "relCoverB_BLTree", 
  "relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
                   "total_NonTree_AbsCover_CONUS")

ggarrange(
# ggarrange(relCover_B_Maps_contemp$relCoverB_totalTree$map, 
#           relCover_B_Maps_contemp$relCoverB_totalTree$map_resids, 
#           relCover_B_Maps_contemp$relCoverB_totalTree$map_future1,
#           relCover_B_Maps_contemp$relCoverB_totalTree$map_future2,
#           relCover_B_Maps_contemp$relCoverB_totalTree$hist,
#           relCover_B_Maps_contemp$relCoverB_totalTree$hist_resids,
#            relCover_B_Maps_contemp$relCoverB_totalTree$hist_future1, 
#           relCover_B_Maps_contemp$relCoverB_totalTree$hist_future2,
#           ncol = 4,
#           nrow = 2,
#           heights = c(1,.35)
#           ),
ggarrange(relCover_B_Maps_contemp$relCoverB_totalHerb$map, 
          relCover_B_Maps_contemp$relCoverB_totalHerb$map_resids, 
          relCover_B_Maps_contemp$relCoverB_totalHerb$map_future1,
          relCover_B_Maps_contemp$relCoverB_totalHerb$map_future2,
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist,
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist_resids,
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_shrub$map, 
          relCover_B_Maps_contemp$relCoverB_shrub$map_resids,
          relCover_B_Maps_contemp$relCoverB_shrub$map_future1,
          relCover_B_Maps_contemp$relCoverB_shrub$map_future2,
          relCover_B_Maps_contemp$relCoverB_shrub$hist,
          relCover_B_Maps_contemp$relCoverB_shrub$hist_resids,
          relCover_B_Maps_contemp$relCoverB_shrub$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_shrub$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_bareGround$map, 
          relCover_B_Maps_contemp$relCoverB_bareGround$map_resids,
          relCover_B_Maps_contemp$relCoverB_bareGround$map_future1,
          relCover_B_Maps_contemp$relCoverB_bareGround$map_future2,
          relCover_B_Maps_contemp$relCoverB_bareGround$hist,
          relCover_B_Maps_contemp$relCoverB_bareGround$hist_resids,
          relCover_B_Maps_contemp$relCoverB_bareGround$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_bareGround$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
# ggarrange(relCover_B_Maps_contemp$relCoverB_NLTree$map, 
#           relCover_B_Maps_contemp$relCoverB_NLTree$map_resids,
#           relCover_B_Maps_contemp$relCoverB_NLTree$map_future1,
#           relCover_B_Maps_contemp$relCoverB_NLTree$map_future2,
#           relCover_B_Maps_contemp$relCoverB_NLTree$hist,
#           relCover_B_Maps_contemp$relCoverB_NLTree$hist_resids,
#           relCover_B_Maps_contemp$relCoverB_NLTree$hist_future1, 
#           relCover_B_Maps_contemp$relCoverB_NLTree$hist_future2,
#           ncol = 4,
#           nrow = 2,
#           heights = c(1,.35)
#           ),
# ggarrange(relCover_B_Maps_contemp$relCoverB_BLTree$map, 
#           relCover_B_Maps_contemp$relCoverB_BLTree$map_resids,
#           relCover_B_Maps_contemp$relCoverB_BLTree$map_future1,
#           relCover_B_Maps_contemp$relCoverB_BLTree$map_future2,
#           relCover_B_Maps_contemp$relCoverB_BLTree$hist,
#           relCover_B_Maps_contemp$relCoverB_BLTree$hist_resids,
#           relCover_B_Maps_contemp$relCoverB_BLTree$hist_future1, 
#           relCover_B_Maps_contemp$relCoverB_BLTree$hist_future2,
#           ncol = 4,
#           nrow = 2,
#           heights = c(1,.35)
#           ),
ggarrange(relCover_B_Maps_contemp$relCoverB_C3Grass$map, 
          relCover_B_Maps_contemp$relCoverB_C3Grass$map_resids,
          relCover_B_Maps_contemp$relCoverB_C3Grass$map_future1,
          relCover_B_Maps_contemp$relCoverB_C3Grass$map_future2,
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist,
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist_resids,
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_C4Grass$map, 
          relCover_B_Maps_contemp$relCoverB_C4Grass$map_resids,
          relCover_B_Maps_contemp$relCoverB_C4Grass$map_future1,
          relCover_B_Maps_contemp$relCoverB_C4Grass$map_future2,
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist,
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist_resids,
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_Forb$map,
          relCover_B_Maps_contemp$relCoverB_Forb$map_resids, 
          relCover_B_Maps_contemp$relCoverB_Forb$map_future1,
          relCover_B_Maps_contemp$relCoverB_Forb$map_future2,
          relCover_B_Maps_contemp$relCoverB_Forb$hist,
          relCover_B_Maps_contemp$relCoverB_Forb$hist_resids,
          relCover_B_Maps_contemp$relCoverB_Forb$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_Forb$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ), 
ggarrange(relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map, 
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_resids,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future1,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future2,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_resids,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future1, 
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ncol = 1
)

Now, we show quantile plots for model predictions of relative cover made with contemporary data. These plots compare final predictions of relative cover to observed relative cover. Quantile plots for relative cover of tree functional types are now show, since they are identical to those for absolute cover.

# First, we need to relativize the predictions and the observations (stored in preds_quantiles)
preds_quantile <- preds_quantile %>% 
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                         absBareGround_CONUS + absShrub_CONUS), # for model predictions
       total_NonTree_AbsCover_obs = (TotalHerbaceousCover + ShrubCover + BareGroundCover) # for observations
       ) %>% 
  # making model predictions relative
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS,
  # making observations relative
         relCoverB_totalTree_Obs = TotalTreeCover,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb_Obs = TotalHerbaceousCover/total_NonTree_AbsCover_obs,
         relCoverB_shrub_Obs = ShrubCover/total_NonTree_AbsCover_obs,
         relCoverB_bareGround_Obs = BareGroundCover/total_NonTree_AbsCover_obs,
         relCoverB_NLTree_Obs = NLTreeCover_abs,
         relCoverB_BLTree_Obs = BLTreeCover_abs,
         relCoverB_C3Grass_Obs = C3GramCover_abs/total_NonTree_AbsCover_obs,
         relCoverB_C4Grass_Obs = C4GramCover_abs/total_NonTree_AbsCover_obs,
         relCoverB_Forb_Obs = ForbCover_abs/total_NonTree_AbsCover_obs,
         )

## exclude trees, since the predictions are the same as the absolute cover shown above
# shrubs
# unique(c(names(coefficients(shrub_CONUS))))
deciles_shrub_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_shrub_Obs", "relCoverB_shrub","prcp",  "prcp_seasonality",  "prcpTempCorr", "sand" ,"coarse",  "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean"  ) %>% drop_na(), 
                                  response_vars = c("relCoverB_shrub_Obs", "relCoverB_shrub"), # name of observations, followed by name of predictions
  pred_vars = c( "prcp",  "prcp_seasonality",  "prcpTempCorr", "sand" ,"coarse",  "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean" ), # for predictors, combine list of predictors for all models 
cut_points = seq(0, 1, 0.01)
                                     )

# bare ground
# unique(c(names(coefficients())))
deciles_bareGround_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_bareGround_Obs", "relCoverB_bareGround", "tmean"  ,"prcpTempCorr"   , "isothermality"   , "annWetDegDays"  ,
"sand"   ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ) %>% drop_na(), 
                                       response_vars = c("relCoverB_bareGround_Obs", "relCoverB_bareGround"), # name of observations, followed by name of predictions
  pred_vars = c( "tmean"  ,"prcpTempCorr"   , "isothermality"   , "annWetDegDays"  ,
"sand"   ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
# total herbaceous cover
#unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS))))
deciles_totalHerb_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_totalHerb_Obs", "relCoverB_totalHerb", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom",   "tmean",  
                "prcp_anom", "clay", "carbon") %>% drop_na(),
                response_vars = c("relCoverB_totalHerb_Obs", "relCoverB_totalHerb"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom",   "tmean",  
                "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

# forbs
# unique(c(names(coefficients(totalTree_F)), names(coefficients(totalTree_GS)), names(coefficients(BLtree_GS)), names(coefficients(BLtree_F))))
deciles_forb_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_Forb_Obs", "relCoverB_Forb", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay", 
                "carbon", "prcp_seasonality",  "annWetDegDays" ,  "prcp_seasonality_anom", "prcpTempCorr_anom" ,  "annWetDegDays_anom") %>% drop_na(), 
                                 response_vars = c("relCoverB_Forb_Obs", "relCoverB_Forb"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay", 
                "carbon", "prcp_seasonality",  "annWetDegDays" ,  "prcp_seasonality_anom", "prcpTempCorr_anom" ,  "annWetDegDays_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

# C3 grass
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C3grass_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_C3Grass_Obs", "relCoverB_C3Grass", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand",  "coarse",  "AWHC" ,  "isothermality_anom",  
                "tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom",  "prcpTempCorr_anom",  "annWetDegDays_anom", "prcp_seasonality" ) %>% drop_na(), 
                                    response_vars = c("relCoverB_C3Grass_Obs", "relCoverB_C3Grass"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand",  "coarse",  "AWHC" ,  "isothermality_anom",  
                "tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom",  "prcpTempCorr_anom",  "annWetDegDays_anom", "prcp_seasonality" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )
# C4 grass 
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C4grass_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_C4Grass_Obs", "relCoverB_C4Grass", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC",  "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na(), 
                                    response_vars = c("relCoverB_C4Grass_Obs", "relCoverB_C4Grass"), # name of observations, followed by name of predictions
  pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC",  "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
                                     )

# shrubs 
quantPlot_shrubs_Obs <- decile_dotplot_pq(df = deciles_shrub_Obs, response= c("relCoverB_shrub_Obs", "relCoverB_shrub"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Shrub Cover: Quantile Plot")

quantPlot_shrubs_Obs <- add_dotplot_inset(quantPlot_shrubs_Obs, df = deciles_shrub_Obs, response = c("relCoverB_shrub_Obs", "relCoverB_shrub"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# bare ground
quantPlot_bareGround_Obs <- decile_dotplot_pq(df = deciles_bareGround_Obs, response= c("relCoverB_bareGround_Obs", "relCoverB_bareGround"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Bare Ground Cover: Quantile Plot")

quantPlot_bareGround_Obs <- add_dotplot_inset(quantPlot_bareGround_Obs, df = deciles_bareGround_Obs, response = c("relCoverB_bareGround_Obs", "relCoverB_bareGround"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# total herbaceous
quantPlot_totHerb_Obs <- decile_dotplot_pq(df = deciles_totalHerb_Obs, response= c("relCoverB_totalHerb_Obs", "relCoverB_totalHerb"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Total Herbaceous Cover: Quantile Plot")

quantPlot_totHerb_Obs <- add_dotplot_inset(quantPlot_totHerb_Obs, df = deciles_totalHerb_Obs, response = c("relCoverB_totalHerb_Obs", "relCoverB_totalHerb"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# forbs
quantPlot_forbs_Obs <- decile_dotplot_pq(df = deciles_forb_Obs, response= c("relCoverB_Forb_Obs", "relCoverB_Forb"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute Forb Cover: Quantile Plot")

quantPlot_forbs_Obs <- add_dotplot_inset(quantPlot_forbs_Obs, df = deciles_forb_Obs, response = c("relCoverB_Forb_Obs", "relCoverB_Forb"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

# C3 grass
quantPlot_c3grass_Obs <- decile_dotplot_pq(df = deciles_C3grass_Obs, response= c("relCoverB_C3Grass_Obs", "relCoverB_C3Grass"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute C3 Grass Cover: Quantile Plot")

quantPlot_c3grass_Obs <- add_dotplot_inset(quantPlot_c3grass_Obs, df = deciles_C3grass_Obs, response = c("relCoverB_C3Grass_Obs", "relCoverB_C3Grass"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

#C4 grass 
quantPlot_c4grass_Obs <- decile_dotplot_pq(df = deciles_C4grass_Obs, response= c("relCoverB_C4Grass_Obs", "relCoverB_C4Grass"), IQR = TRUE,
                        CI = FALSE
                        ) + ggtitle("Absolute C4 Grass Cover: Quantile Plot")

quantPlot_c4grass_Obs <- add_dotplot_inset(quantPlot_c4grass_Obs, df = deciles_C4grass_Obs, response = c("relCoverB_C4Grass_Obs", "relCoverB_C4Grass"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)

## now, display the figures
quantPlot_totHerb_Obs

quantPlot_shrubs_Obs

quantPlot_bareGround_Obs

quantPlot_c3grass_Obs

quantPlot_c4grass_Obs

quantPlot_forbs_Obs

Below, we show the RMSE and bias of the final estimates of relative cover for each functional type (except trees) in comparison to the observed data.

long_absPreds <- preds_quantile %>% 
  mutate(uniqueID = 1:nrow(preds_quantile)) %>% 
  select(uniqueID, relCoverB_totalHerb, relCoverB_shrub,             
relCoverB_bareGround, relCoverB_C3Grass,           
relCoverB_C4Grass, relCoverB_Forb) %>% 
  rename(TotalHerb = relCoverB_totalHerb, 
         Shrub = relCoverB_shrub, 
         BareGround = relCoverB_bareGround, 
         C3gram = relCoverB_C3Grass, 
         C4gram = relCoverB_C4Grass, 
         Forb = relCoverB_Forb) %>% 
  pivot_longer(cols = c(TotalHerb, Shrub, BareGround, C3gram, C4gram, Forb), 
                names_to = "names", 
               values_to = "preds_values"
               ) %>% 
  select(uniqueID, names, preds_values)

long_absObs <- preds_quantile %>% 
  mutate(uniqueID = 1:nrow(preds_quantile)) %>% 
  select(uniqueID, relCoverB_totalHerb_Obs, relCoverB_shrub_Obs,  relCoverB_bareGround_Obs, relCoverB_C3Grass_Obs, relCoverB_C4Grass_Obs,  relCoverB_Forb_Obs) %>% 
  rename(TotalHerb = relCoverB_totalHerb_Obs, 
         Shrub = relCoverB_shrub_Obs, 
         BareGround = relCoverB_bareGround_Obs, 
         C3gram = relCoverB_C3Grass_Obs, 
         C4gram = relCoverB_C4Grass_Obs, 
         Forb = relCoverB_Forb_Obs) %>% 
  pivot_longer(cols = c(TotalHerb, Shrub, BareGround, C3gram, C4gram, Forb), 
                names_to = "names", 
               values_to = "obs_values"
               ) %>% 
  select(uniqueID, names, obs_values)

longObsPreds <- long_absObs %>% 
  left_join(long_absPreds)

## calculate bias
(summaryTableAbs <- longObsPreds %>% 
  mutate(resids = obs_values - preds_values) %>% 
  group_by(names) %>% 
  dplyr::summarise(bias = mean(resids, na.rm = TRUE),
            RMSE = round(sqrt(mean(resids^2, na.rm = TRUE)),4)
            ) %>% 
  slice(8, 9, 7, 2, 6, 1, 3, 4, 5) %>% kable(
  format = "html", 
  col.names = c("Model Name", "Model bias mean(obs. - pred.)", "model RMSE")
))

Model Name	Model bias mean(obs. - pred.)	model RMSE
C3gram	0.0057089	0.1750
TotalHerb	-0.0022415	0.1421
BareGround	0.0044805	0.1227
C4gram	-0.0005822	0.1017
Forb	-0.0185035	0.1360
Shrub	-0.0022389	0.1043

Finally, we show how these model predictions made using contemporary climate and weather data change over time from 2010 to 2020.

# calculate relativized cover (version w/ all functional groups summing to 1)
predsOverTime <- predsOverTime %>% 
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                         absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

# make into a long format
predsOverTime_long <- predsOverTime %>% 
  select(year,x, y, spatialID, uniqueID, total_NonTree_AbsCover_CONUS:relCoverB_Forb) %>% 
  pivot_longer(cols = c(total_NonTree_AbsCover_CONUS:relCoverB_Forb)) %>% 
  mutate(year = as.integer(year))

# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>% 
  group_by(uniqueID, name) %>% 
  summarize(cv = sd(value)/mean(value))
## add x,y 
predsOverTime_cv <- predsOverTime_cv %>% 
  left_join(predsOverTime %>% select(uniqueID, x, y))

# visualize change over time for a selection of points
coverOverTimeForSelectLocations <- predsOverTime_long %>% 
  filter(uniqueID %in% 1:50) %>% 
  filter(name %in% c( #"total_NonTree_AbsCover_CONUS",  "relCoverB_totalTree",    
#"relCoverB_totalHerb", 
    "relCoverB_shrub"    , "relCoverB_bareGround" , "relCoverB_NLTree",  "relCoverB_BLTree",       
"relCoverB_C3Grass" , "relCoverB_C4Grass",   "relCoverB_Forb" )) %>% 
  ggplot() + 
  facet_wrap(~uniqueID) + 
  geom_area(aes(x = year, y = value, fill = name)) + 
  geom_line(data = predsOverTime_long %>% 
              filter(uniqueID %in% 1:50) %>% 
              filter(name %in% c("relCoverB_totalHerb", "relCoverB_totalTree")),
            aes(x = year, y = value, linetype = name), col = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90)) + 
  labs(title = "Change in modeled relative cover (all Non-tree types sum to 1) by functional type from \n 2010 to 2020 for a subset of locations")
mapOfSelectLocations <- predsOverTime_long %>% 
    select(x, y, uniqueID) %>% 
    filter(uniqueID %in% 1:50) %>% 
    unique() %>% 
    ggplot() + 
    geom_sf(data = cropped_states, aes()) + 
    geom_point(aes(x, y)) + 
    geom_label_repel(aes(x, y, label = uniqueID)) + 
    labs(title = "Locations in above plot")
    

ggarrange(coverOverTimeForSelectLocations, mapOfSelectLocations, ncol = 1, heights = c(1,.5))

predsOverTime_cv %>%  
  ggplot() + 
  geom_boxplot(aes(x = name, y = cv, col = name)) +
  theme_minimal() + 
  labs(title = "Coefficients of variation of model-predicted  relative cover (all Non-tree types sum to 1) \n at a given location for each functional group from 2010-2020",
        x = "functional group",
       y = "coefficient of variation") +
  theme(axis.text.x = element_text(angle = 90))

# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>% 
 # rename(geometry = Shape)
 
ggplot() + 
  facet_wrap(~name) + 
  geom_sf(data = cropped_states, aes()) + 
  stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
  geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) + 
  guides(fill = guide_legend(title = "CV")) + 
  theme_minimal() + 
  labs(title = "The spatial distribution of the coefficients of variation of \nmodel-predicted cover from 2010-2020 for 1000 random locations",
       subtitle = "Note that values are aggregated accross space to enhance readability")